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\  /  ^  ABSTRACT 

^his  report  evaluates  superresolution  direction  finding  methods  for  land  tactical 
narrowband  \4JFy-WiF  purposes.  The  DF  methods  are  described  and  evaluated  in  depth 
using  a  common  theoretical  framework  beginning  with  classical  methods  and  ending  with 
modern  day  eigenanalysis  based  methods.  Based  on  this  analysis,  superresolution  methods 
can  be  described  in  terms  of  a  five  step  procedure  where  each  step  has  a  unique  purpose  in 
the  estimation  process.  Different  methods  vary  according  to  the  operations  performed  at 
each  step.  This  is  useful  in  analyzing  and  comparing  the  theoretical  performance  of  various 
DF  methods  under  the  same  conditions.  The  results  from  simulations  are  also  included  to 
support  the  theoretical  evaluations.  i- '  -.u 


RESUME 

Ce  rapport  contient  une  evaluation  des  techniques  de  radiogoniometrie  a 
superresolution  pour  des  signaux  tactiques  terrestres  VHF/UHF  a  bande  passante  etroite. 
Les  methodes  de  radiogoniometrie  sont  decrites  et  evaluees  en  detail  a  partir  d’une  theorie 
unifiee  valide  pour  les  methodes  classiques  et  les  methodes  modernes  utilisant  des  valeurs 
propres.  II  est  possible,  a  partir  de  cette  theorie,  de  decrire  les  methodes  de  superresolution 
a  I’aide  d’une  procedure  en  cinq  etapes  ou  chaque  etape  a  un  but  unique  dans  le  processus 
d’evaluation.  Les  differentes  methodes  se  differencient  par  les  operations  effectuees  a 
chaque  etape.  II  s’agit  d’une  caracteristique  utile  pour  I’analyse  et  la  comparaison,  dans  les 
memes  conditions,  des  performances  theoriques  de  diverses  methodes  de  radiogoniometrie. 
Les  evaluations  theoriques  presentees  sont  etayees  par  des  simulations. 
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EXECUTIVE  SUMMARY 


A  considerable  amount  of  work  on  superresolution  radio  direction  finding  (DF) 
methods  has  been  reported  in  open  literature  sources  over  the  last  twenty  years.  However, 
due  to  the  variability  in  approaches,  it  is  difficult  to  make  definitive  statements  as  to  the 
relative  performance  and  merits  of  the  various  methods  described.  In  this  report  a  common 
theoretical  framework  is  used  to  describe  the  various  methods.  The  theoretical  discussion 
presented  in  this  report  is  based  on  the  analysis  as  firjt  presented  in  the  open  literature  as 
well  as  new  analysis.  The  new  analysis  was  required  to  ensure  a  unified  approach  in  the 
theoretical  discussion  as  well  as  provide  a  firm  theoretical  basis  for  techniques  which  have 
been  taken  for  granted  in  the  open  literature,  but  rarely  discussed. 

Simulations  are  also  presented  to  illustrate  relative  performance  of  the  various  DF 
methods  using  the  same  input  data.  Since  this  report  was  written  in  support  of  tactical 
VHF/UHF  DF,  the  simulations  deal  predominantly  with  multipath  (fully  coherent)  signals 
since  multipath  has  been  perceived  to  be  the  most  serious  source  of  error  in  current 
operational  DF  systems. 

To  keep  this  report  to  a  manageable  size,  only  representative  DF  methods  based 
on  classical  spectral  estimation  techniques  and  those  which  have  evolved  from  these 
techniques  have  been  described.  For  example,  this  includes  the  Bartlett  method.  Linear 
Prediction,  MUSIC,  and  other  similar  estimators,  but  does  not  include  the  phase 
interferometer.  Esprit  and  Maximum  Likelihood  methods.  Additionally,  analysis  has  been 
restricted  to  uniform  linear  arrays  and  white  Gaussian  noise  statistics.  Despite  these 
limitations,  the  discussion  in  this  report  is  reflective  of  the  main  body  of  research  in  DF 
over  the  last  twenty  years  and  before  this  report  was  written. 

DF  methods  can  be  described  in  terms  of  three  filter  models,  namely,  the  movin'^ 
average  filter,  the  autoregressive  filter,  and  the  autoregressive  moving  average  filter.  In 
terms  of  accuracy  and  computational  speed,  methods  based  on  the  all  pole  or  autoregressive 
filter  model  have  shown  the  greatest  promise  and  for  this  reason  were  the  main  focus  of  this 
report. 


It  is  shown  that  the  autoregressive  filter  based  methods  can  be  generalized  as  a 
five  step  procedure.  These  steps  include; 

1.  Estimation  of  the  autocorrelation  matrix. 

2.  Division  of  the  autocorrelation  matrix  into  a  signal  and  noise  subspace. 

3.  Generation  of  an  enhanced  inverse  autocorrelation  matrix. 

4.  Estimation  of  the  all  pole  filter  coefficients. 

5.  Estimation  of  the  signal  bearings  from  the  filter  coefficients. 

Each  of  these  steps  is  the  subject  of  separate  sections  in  this  report 

In  terms  of  the  resultant  DF  accuracy,  and  given  the  assumptions  and  restrictions 
made  in  this  report,  the  root-MUSIC  and  root-Minimum  Norm  methods  represent  two  of 
the  best  methods  for  DF.  These  methods  are  far  superior  to  classical  techniques  but  are  in 
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no  way  optimum  -  more  research  is  still  required  in  each  of  the  steps  listed  previously. 
Additionally,  in  practice,  many  of  the  assumptions  and  restrictions  used  in  this  report  are 
often  untrue,  so  that  other  areas  of  research  include: 

1.  Extension  of  techniques  used  for  uniform  linear  arrays  to  arbitrary  array 
geometries,  or  development  of  comparable  techniques. 

2.  Sensor  calibration. 

3.  Estimation  of  noise  statistics  in  an  unknown  coloured  noise  environment 

4.  Modifications  to  the  signal  model  for  actual  multipath  environments  and 
mixed  coherent /noncoherent  signal  environments. 

5.  Faster  algorithms  for  realtime  implementation. 

Solutions  to  the  above  problems  will  be  required  before  the  promises  of  superresolution 
methods  can  be  fully  realized. 
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1.0  INTRODUCTION 


A  considerable  amount  of  research  and  development  has  been  carried  out  over  the 
last  twenty  years  in  advanced  spectral  estimation  methods  which  has  led  directly  to  the 
development  of  advanced  radio  direction  finding  (DF)  methods.  The  methods  of  interest 
are  those  which  are  intended  to  operate  against  conventional  narrowband  VHF/UHF  radio 
signals  (i.e.,  CW,  AM,  FM,  and  SSB  with  bandwidths  less  than  200  kHz)  in  a  land-based 
tactical  environment.  This  report  discusses  these  methods  and  their  potential  for 
improving  the  capabilities  of  current  conventional  ^  tical  DF  systems  (discussed  in 
reference  [1-1]). 

Each  of  the  methods  are  discussed  in  terms  of  their  development  philosophy  and 
their  ability  to  deal  with  various  error  sources  including  noise,  co-channel  interference,  and 
multipath.  Although  the  subject  of  advanced  DF  is  not  limited  to  theory,  a  fundamental 
understanding  of  the  capabilities  of  these  methods  is  essential  if  a  successful  operational 
system  is  to  be  built. 

1.1  SYSTEM  OVERVIEW 

For  the  sake  of  simplicity,  the  sensor  systems  discussed  in  this  report  are  assumed 
to  consist  of  a  linear  array  of  N  vertical  dipole  antennas,  all  of  which  are  connected  to  N 
separate  perfectly  matched  (in  both  gain  and  phase)  and  like-numbered  receiver  channels  as 
shown  in  Figure  1.1. 


FIGURE  1.1:  Block  diagram  of  a  multi-channel  DF  system 


In  the  receiver,  each  channel  is  filtered  to  the  correct  bandwidth,  mixed  using  a 
common  reference  signal  (at  the  filter  center  frequency)  to  generate  an  in-phase  and 
quadrature  baseband  representation  of  the  antenna  signal,  which  are  then  digitized.  The 
sampling  rate  is  assumed  to  be  twice  the  bandwidth  (i.e.  the  Nyquist  rate).  The 
bandwidth  of  the  signal  is  assumed  to  be  sufficiently  narrow  so  that  the  modulation 
envelope  does  not  change  during  the  time  taken  for  a  signal  to  traverse  the  array  when  the 
signal  is  incident  from  boresight. 


1 


The  receiver  outputs  are  used  as  the  inputs  to  the  bearing  processor  which  are 
represented  by  the  complex  valued  parameters  lo,  2:2,  xv-i  where  the  subscript 
represents  the  channel  number.  The  bearing  processor  inputs  are  assumed  to  be 
undistorted  baseband  representations  of  the  antenna  signals.  In  this  context  the  values 
lo,  xi,  2:2,  ...,  i.v-i  are  referred  to  as  the  sensor  data.  A  single  sample  (alternately  called  a 
"snapshot"  in  the  open  literature)  is  defined  as  the  set  of  values  lo,  x\,  12,  2:,v-i  measured 

at  a  specific  instance  in  time.  To  distinguish  sensor  values  from  the  same  channel  n  but 
measured  at  different  points  in  time  (i.e.  different  time  samples),  an  alternate 
representation  is  used,  namely,  Xn{t),  where  t  represents  the  sample  index. 

Throughout  this  report,  the  parameters  M,  iV,  and  Tare  used  exclusively  to 
represent  the  number  of  signals,  the  number  of  sensors,  and  the  number  of  sensor  data 
samples,  respectively.  Other  definitions  and  conventions  followed  throughout  this  report 
are  described  in  the  glossary. 

1.2  SPFICTRAL  ESTIMATION  AND  RADIO  DIRECTION  FINDING 

The  application  of  spectral  estimation  concepts  to  radio  direction  finding  is 
illustrated  using  Figure  1.2.  In  this  representation  a  single  radio  signal  in  a  noiseless 
environment  traverses  an  antenna  array  with  a  bearing  in  azimuth  given  by  <p.  The 
assumption  is  made  here  (and  throughout  the  rest  of  the  report)  that  the  transmitter  is  far 
enough  from  the  antenna  array  so  ■‘hat  the  wavefront  is  planar.  Additionally,  since  at 
VFIF/UHF  signal  propagation  is  generally  limited  to  ground  waves,  elevation  angles  are 
always  considered  to  be  0  degrees. 


FIGURE  I  2:  Radio  direction  finding  sensor  array 


The  amplitude  and  phase  measured  at  sensor  n  in  Figure  1.2  is  represented  at 
baseband  by  the  complex  value, 


(1.1) 


where  c{t)  is  the  complex  baseband  modulation  envelope  measured  at  sensor  0,  d  is  the 
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sensor  spacing,  and  u  is  the  spatial  frequency  defined  by 

6J  =  ^cos((^)  rad/m.  (1.2) 

Here  A  represents  the  signal  wavelength.  The  term  spatial  frequency  is  used  to  denote  that 
the  frequency  is  measured  with  respect  to  position  rather  than  vnth  respect  to  time.  For 
more  complex  antenna  configurations,  the  spatial  frequency  may  still  be  calculated  using 
equation  (1.2)  by  designating  an  arbitrary  baseline  as  the  reference  baseline  and  then 
making  all  angular  measurements  with  respect  to  it. 

If  a  number  of  signals  of  different  bearings  impinge  upon  the  sensor  system  in  a 
noisy  environment,  then  the  spatial  signal  x„  is  the  sum  of  the  individual  signals,  that  is, 

Xn  =  (1.3) 

m=  I 

where  Af  represents  the  total  number  of  signals  and  r/„(t)  is  the  noise  measured  at  sensor  n. 
Spectral  estimation  methods  can  be  used  to  resolve  the  various  components  in  equation 
(1.3)  by  first  generating  the  spatial  power  spectral  density  function.  By  taking  advantage 
of  the  relationship  between  spatial  frequency  and  bearing,  this  spectrum  is  converted  to  a 
DF  spectrum  which  gives  the  power  density  versus  bearing.  The  location  of  the  peaks  (local 
maximums)  in  the  DF  spectrum  are  then  used  as  estimates  for  the  signal  bearings. 

The  spectral  estimation  approach  has  led  to  the  development  of  a  number 
sophisticated  direction  finding  methods  whose  origins  are  based  on  research  done  in  fields 
as  diverse  as  geophysics,  digital  communications,  underwater  acoustics,  and  speech 
analysis.  In  many  of  these  applications  spectral  estimation  methods  are  used  to  process 
time  series  (or  temporal)  data,  whereas  in  the  direction  finding  case  these  methods  are  used 
to  process  position  series  (or  spatial)  data. 

Mathematically  there  are  some  important  differences  between  temporal  frequency 
estimation  methods  and  spatial  frequency  estimation  methods.  One  difference  is  that 
position  is  a  three  dimensional  parameter  while  time  is  one  dimensional.  Unless  otherwise 
indicated,  the  DF  methods  in  this  report  are  discussed  in  terms  of  a  single  antenna  baseline 
with  uniform  antenna  spacing.  In  this  special  case,  position  is  restricted  to  one  dimension 
and  the  correspondence  between  the  temporal  and  spatial  methods  is  one  to  one.  In 
general,  results  for  single  baselines  can  be  extended  to  multiple  baselines,  if  care  is  used. 

Another  important  difference  between  temporal  and  spatial  methods  is  that  in 
spatial  frequency  estimation  the  time  dimension  has  no  equivalence  in  temporal  frequency 
estimation.  This  extra  dimension  provides  information  which  is  useful  when  dealing  with 
uncorrelated  signals  or  temporal  noise.  It  is  also  easier  and  a  lot  less  expensive  to  sample 
with  respect  to  time  than  with  respect  to  position. 

The  layout  of  this  report,  where  applicable,  follows  the  historical  development  of 
the  DF  methods  discussed,  starting  with  classical  methods  (Section  2),  followed  by  model 
based  methods  (Sections  4-5),  then  on  to  the  more  advanced  eigenanaJysis  methods 
(Sections  6-8),  and  finally  ening  with  a  discussion  of  the  rooting  method  (Section  9). 
Additional  sections  have  also  been  included  for  the  discussion  of  methods  which  can  be 
considered  separately  from  DF,  but  are  nonetheless  critical  to  the  success  of  advanced  DF 
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methods.  This  includes  autocorrelation  matrix  estimation  (Section  3),  and  model  order 
determination  (Section  9). 

The  methods  discussed  in  this  report  are,  for  the  most  part,  based  on  the 
autocorrelation  matrix  formulations  since  this  simplifies  comparisons.  In  many  cases 
variants  of  these  methods  exist  which  are  computationally  more  attractive,  but  are  not 
discussed  since  they  are  theoretically  identical. 
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2.0  THE  CLASSICAL  APPROACH 


The  exact  DF  spectrum  can  be  calculated  by  computing  the  spatial 
autocorrelation  sequence  and  then  taking  the  Fourier  transform  of  this  sequence.  This  is 
called  the  indirect  method.  Alternatively,  the  exact  DF  spectrum  can  be  calculated  by 
taking  the  Fourier  transform  of  the  input  data  sequence  and  squaring  the  absolute  values  of 
the  transformed  sequence.  This  is  called  the  direct  method. 

In  either  the  direct  or  indirect  methods,  the  assumption  is  made  that  there  is  an 
infinite  amount  of  data  available.  Under  realistic  conditions,  this  clearly  will  not  be  the 
case,  so  instead  approximations  must  be  used  which  result  in  DF  spectrums  which  are  only 
estimations  of  the  true  function.  Additionally,  under  these  conditions,  the  two  methods 
may  produce  different  results. 

2.1  THE  INDIRECT  METHOD  OF  ESTIMATING  THE  DF  SPECTRUM 

The  DF  spectrum  can  be  defined  as  the  Fourier  transform  of  the  spatial 
autocorrelation  sequence  and  is  given  by, 


S{<P)  =  X  (2.1) 

m--  <30 

where  the  spatial  frequency  ui  is  related  to  the  bearing  angle  (f>  by  equation  (1.2),  and  d 
is  the  spacing  between  consecutive  antennas.  The  definition  and  estimation  of  the 
autocorrelation  sequence  r-i(m)  is  discussed  in  the  next  two  sections. 

2.1.1  The  Autocorrelation  Sequence 

The  statistical  definition  of  the  autocorrelation  of  a  spatially  sampled  random 
process,  i*,  at  two  different  indices  m  and  n  is  defined  here  as 

r„(m,7i)  =  I*}  (2.2) 

where  E{y}  is  the  mean  or  expected  value  of  the  process  represented  by  y.  If  the 
sampled  process  is  wide  sense  stationary,  then  its  mean  is  constant  for  ail  indices,  or 

E{xm}  =  E{xr,},  (2.3) 

and  its  autocorrelation  will  depend  only  on  the  difference  n-m,  or 


E{x,n  ij}  =  E{x,„,k  In**}.  (2-4) 

where  k  is  any  arbitrary  integer.  Under  these  conditions,  an  autocorrelation  sequence 
may  be  defined  where  each  element  is  given  by 

r^j{m)  =  E{xr,,mxt}-  (2.5) 

The  index  value  m  is  called  the  autocorrelation  lag. 

In  direction  finding,  the  wide  sense  stationary  condition  means  that  all  signal 
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sources  must  be  located  in  the  farfield  of  the  antenna  array  (i.e.  the  wavefront  arriving  at 
the  anteima  array  produced  by  any  one  source  is  planar)  and  the  first  and  second  order 
statistics  of  the  noise  (mean  and  variance)  do  not  change  with  respect  to  position  and  time. 

2.1.2  The  Unbiased  Estimator 

A  more  practical  definition  of  the  exact  autocorrelation  sequence  than  the  one 
defined  by  equation  (2.5),  which  takes  advantage  of  the  fact  that  the  input  data  is  sampled 
in  both  position  and  time,  is  given  by 


r  s 

r„{Tn)  -  ^Hm^  2T+T  2  277+T  2  (2-6) 

’  «=-T 

where  T  represents  the  number  of  time  samples  and  N  represents  the  number  of  sensors 
(these  definitions  for  T  and  N  are  used  throughout  the  rest  of  this  report).  Here  t  represents 
the  temporal  index  and  n  the  spatial  index.  The  Fourier  transform  of  this  sequence  is  the 
exact  spectrum.  Unfortunately,  under  realistic  conditions,  the  number  of  data  values  is 
finite  and  consequently  the  value  of  Nis  also  finite,  rather  than  infinite.  As  a  result,  it  will 
only  be  possible  to  estimate  the  autocorrelation  sequence. 

One  method  of  estimating  the  sequence  is  to  modify  equation  (2.6)  to  give, 

r-i  y-m-i 

r„(m)  =-f^  Tra  2  for  0  <  m  <  iV,  (2.7) 

1=0  n  =  0 


=  0  for  m  >  A,  (2.8) 

and, 


rxJ[m)  =  rt^-m)  for  m  <  0.  (2.9) 

Equation  (2.7)  is  an  unbiased  estimate  of  the  autocorrelation  sequence,  since  the  expected 
vadue  of  the  estimated  autocorrelation  lag  is  the  true  value. 

Since  only  estimates  of  the  autocorrelation  sequence  are  available,  the  Fourier 
transform  of  this  sequence  will  be  an  estimate  of  the  DF  spectrum.  That  is. 


5(0)  =  ^  (2.10) 

m»l  -W 

where  the  values  of  r„{m)  are  computed  using  equations  (2.7)  and  (2.9). 

2.1.3  The  Biased  Estimator 

One  difficulty  with  generating  the  autocorrelation  sequence  using  equation  (2.7), 
especially  if  only  a  few  samples  are  available,  is  that  the  variances  of  the  estimated 
autocorrelation  sequence  values  increase  as  the  absolute  lag  index  increases,  since  there  are 
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less  and  less  data  values  available  for  averaging.  As  a  consequence  of  this,  it  is  possible  for 
the  estimated  sequence  to  be  non-realizable.  For  example,  the  estimate  of  any  value  in  the 
autocorrelation  sequence  (except  at  lag  0)  may  result  in  a  value  larger  than  the  one  at  lag 
0,  which  is  not  physically  possible. 

To  overcome  this  problem,  a  triangular  weighting  function  (window)  is  applied 
which  results  in  the  modified  estimate  given  by 

T-  1  H-m-l 

r„(m)  =  ^  Z 

t-0  n=0 

This  is  a  biased  estimate  of  the  autocorrelation  sequence,  since  the  expected  value  of 
estimated  autocorrelation  lag  is  not  the  true  value,  except  for  lag  0.  However,  this  estimate 
always  results  in  a  physically  realizable  autocorrelation  sequence.  Additionally,  as  the 
number  of  values  N  approaches  infinity,  the  estimated  autocorrelation  sequence 
approaches  the  true  sequence.  In  this  context,  the  estimator  is  asymptotically  unbiased. 

A  second  advantage  of  the  biased  autocorrelation  sequence  estimate,  compared  to 
the  unbiased  estimate,  is  a  lower  variance  since  the  effect  of  the  larger  lag  elements,  which 
are  estimated  from  fewer  product  terms  and  therefore  have  greater  variance,  are 
suppressed.  The  disadvantage  is  that  windowing  degrades  the  resolution  of  the 
corresponding  DF  spectrum  estimate  (which  is  computed  using  equations  (2.10)  and  (2.11)) 

2.2  THE  DIRECT  METHOD  OF  ESTIMATING  THE  DF  SPECTRUM 

An  alternate  method  (called  the  direct  method  here)  of  computing  the  DF 
spectrum  is  to  square  the  absolute  magnitude  of  the  Fourier  transform  calculated  directly 
from  the  data,  that  is, 

(2.12) 

where  X{4))  represents  the  spatial  Fourier  transform  of  the  data. 

If  only  a  finite  amount  of  data  is  available,  X{(p)  can  be  approximated  by 

m=  0 

where  d  is  the  antenna  spacing.  Using  this  approximation,  the  calculated  value  of  5(^) 
will  be  an  estimate  of  the  true  value.  Substituting  equation  (2.13)  into  (2.12)  gives 


k<l>)  =  ^ 

m=0  m*0 

The  spectrum  generated  for  a  random  process  using  ^nation  (2.14)  is  statistically 
unstable  [2-1].  That  is,  as  N  increases,  the  variance  of  the  estimated  vmue  of  S(0)  does  not 
approach  zero,  but  rather  becomes  proportional  to  the  square  of  the  true  value  of  S{((>). 

The  method  used  to  get  around  this  difficulty  is  to  compute  the  estimated  spectrum  for 
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several  different  time  samples,  then  average  the  spectrums  together  to  give  an  average 
version.  That  is, 


k<l>)  =  7^  2  ^  .  (2.15) 

t=0  m=0  m=0 

2.3  MATRIX  REPRESENTATION  OF  THE  CLASSICAL  ESTIMATORS 

2.3.1  The  Direct  Method 

Inspection  of  equation  (2.13),  the  Fourier  transform  of  the  sensor  data,  suggests  a 
compact  matrix  representation  which  is, 


XW  = 


(2.16) 


where  e  and  Xj  are  both  N  element  column  vectors.  The  vector  Xi  is  the  sensor  data  vector 
measured  at  a  time  t,  and  is  defined  as 


a^(<)  ■ 

Xi{t) 

xit) 


(2.17) 


The  N  element  vector  e  is  referred  to  as  the  steering  vector  and  is  given  by 

1 

*jud 

e 

e  =  e 

e 


(2.18) 


The  estimated  DF  spectrum,  using  the  direct  method  (equation  (2.12)),  is  then 


(2.19) 


Following  the  example  of  equation  (2.13)  and  incorporating  time  averaging  to  generate  a 
more  stable  estimate  gives, 


•S'(<^)  =  7^  ^  e“x<X(fe, 


t  =0 


(2.20) 


or. 
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(2.21) 


k<li)  = 

where 

a;o(0),  2^(1),  io(2),  ...,  10(7-1) 
a:i(0)»  *1(1)/  *i(2), xi(r-l) 

X  =  12(0),  12(1),  12(2),...,  12(7-1)  .  (2.22) 


L  iv-i(O),  iv-i(l),  Xs.i{2),  iv-i(7-l)  J 

Equation  (2.21)  is  called  the  Bartlett  estimator  [2-2]  (note  that  since  for  bearing  estimation 
purposes  only  the  shape  of  the  spectrum  is  important  the  factor  1/lVcan  be  dropped  from 
this  expression). 

2.3.2  The  Indirect  Method 

An  alternate  representation  of  the  Bartlett  estimator  given  by  equation  (2.21)  can 
be  formed  be  letting  the  matrix  product  XX^  be  represented  by  a  single  matrix  so  that, 

=  (2.23) 

The  elements  of  the  matrix  R  are  given  by 

*;(0-  (2-24) 

( >0 

The  maximum  stability  in  the  estimate  of  the  spectrum  occurs  as  the  number  of 
samples,  7,  is  increased  to  infinity.  In  this  limiting  case  the  elements  of  the  matrix 
become, 

E{rij}  =  E{a;  x*}  =  r^j{i-j),  (2.25) 

which  are  the  values  of  the  spatial  autocorrelation  sequence.  As  a  result,  R  is  called  the 
autocorrelation  matrix  and  can  be  defined  in  terms  of  the  estimated  autocorrelation  matrix 
as, 

R  =  E{R}.  (2.26) 

The  matrix  R  is  also  called  the  covariance  matrix  since  Xij  also  represents  the  covariance 
between  two  time  varying  processes  defined  by  Xi{t)  and  z/t). 

If  the  matrix  calculations  are  performed  for  equation  (2.23)  the  result  can  be 
expressed  as. 
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(2.27) 


Af-  1  «-  1 


5(0)  -  ^  X  Z 

ms 0  n * 0 


-jumd  “  *jund 
^  ‘^mn 


Replacing  r^i  with  equation  (2.24)  gives 

H-l  N-l  T- 1 


i^^)  =  w  Z  2  T  ^ 


♦  n-mjd 


m=  0  n  =  0  <  =0 


(2.28) 


Grouping  terms  with  the  same  Fourier  coefficient  and  expanding, 

^0)  =  7^[(io(f)  ili(f))  xlit)  +  x,{t) 

+  (l0(t)  lAf-3(t)  +  Zl(t)  lAf-2(^)  +  2:2(0  ZAr-l(0)e*^‘^^*^‘* 

+  ...  +(z;v-i(0  io(0)  (2.29) 

and  rewriting, 

m=l-y  <*0  nsO 

Inspection  of  the  quantity  in  the  square  brackets  above  reveals  it  to  be  the  biased 
autocorrelation  sequence  estimator  defined  earlier  by  equation  (2.11).  From  this,  it  is 
apparent  that  the  use  of  the  autocorrelation  matrix  in  equation  (2.23)  is  equivalent  to  the 
direct  method  if,  and  only  if,  the  biased  autocorrelation  sequence  estimator  is  used. 

2.4  DATA  WINDOWING 


There  are  several  consequences  of  the  spatially  limited  data  available.  For  one, 
the  spatial  frequency  resolution  of  this  technique  is  limited  to  roughly  the  reciprocal  of  the 
antenna  baseline,  or  worse.  Another  consequence  is  that  the  sudden  transition  to  zero  of 
the  unknown  values  in  the  data  sequence  has  the  effect  of  causing  large  side  lobes  (Gibbs 
phenomenon)  in  the  spatial  frequency  spectrum.  Window  functions  can  be  applied  which 
taper  the  data  at  the  start  and  end  of  the  data  sequence  to  reduce  the  sidelobe  problem.  As 
a  consequence,  the  data  is  modified  to  become 


Vn  —  ®n  (2.31) 

where  a„  is  the  weighting  coefficient.  The  modified  value  is  used  in  place  of  Xn  in 
any  of  the  estimation  tec^ques  discussed  so  far. 

In  terms  of  the  matrix  equations  discussed  previously,  the  modification  given  in 
equation  (2.31)  results  in  a  new  data  matrix  Y  (which  replaces  X)  of  the  form. 
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Y  = 


rr 


aoa:o(0),  aoXo(l),  aoio(2),  aoXo(r-l) 

aiii(O),  0111(1),  aiii(2),  oiii(r-l) 

02X2(0),  02X2(1),  02X2(2),  02X2(!r-l) 

[  Ow.iljv.i(0),  OiV-lX;v-i(l),  OjVLiIjv-l(2),  Ojv-ilAr-i(r-l)  J 


(2.32) 


The  effects  of  using  different  types  of  data  windows  are  given  in  reference  [2-3]. 
The  application  of  the  weighting  function  has  a  price,  however,  as  sidelobe  reduction  is 
achieved  at  the  expense  of  resolution. 


2.5  EXTENSION  TO  NONLINEAR  AND  NONUNIFORM  ARRAYS 


Inspection  of  the  steering  vector,  e,  reveals  that  it  contains  the  (conjugate) 
spatial  Fourier  transform  coefficients  corresponding  to  each  antenna  in  the  array.  In 
spatial  terms,  the  steering  vector  represents  a  sinusoidal  signal  propagating  through  space 
in  a  known  direction.  Each  element  in  the  vector  defines  the  amplitude  and  phase 
parameters  of  the  signal  at  each  antenna  position  relative  to  antenna  position  0. 

For  a  sinusoidal  wave  propagating  in  a  known  direction,  the  steering  vector 
elements  can  be  defined  for  any  antenna  position  as 


‘w  (x„cos(^)cos(i)  ♦  y„sin(^)cos(lli)  ♦  z„sin(tl’)) 

e„  =  e  (2.33) 

where  the  antenna  position  is  specified  in  terms  of  the  cartesian  coordinates  (u,  yn,  Zn), 
antenna  0  is  located  at  the  origin,  the  plane  is  the  surface  of  a  flat  Earth,  ana  en  is  the 
element  of  the  steering  vector.  Additionally  the  elevation  bearing  angle  ^  has  also  been 
added  for  the  more  gener^  case  where  the  signal  elevation  angle  must  also  be  determined 
(this  changes  the  bearing  estimation  procedure  from  a  one  dimensional  search  in  azimuth 
to  a  two  dimensional  search  in  both  azimuth  and  elevation). 

Equation  (2.33)  may  be  written  in  a  more  compact  form,  namely. 


e„=  e 


2ir  T 


(2.34) 


where  Pn  is  the  vector  representing  the  position  of  the  nf*  antenna  and  is  defined  in 
Cartesian  coordinates  as 


P«  = 


yn  , 

2n. 


and  B  is  the  signal  direction  vector  defined  by 


(2.35) 
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'  cos(^)cos(^) 
B=  sin(0)sin(^) 
sin(^) 


(2.36) 


For  a  uniform  linear  array,  where  the  array  baseline  is  aligned  with  the  x  axis 
(i.e.  P„  =  [nd,  0,  0]t)  equation  (2.34)  reduces  to  a  more  familiar  form  (see  equation  (2.18)) 
where, 


(2.37) 


2.6  BEAMFORMING 

The  directional  response  of  the  spectral  estimators  described  to  this  point  are 
equivalent  to  beamforming  [2-4].  For  example,  the  general  form  for  the  output  of  a 
beamformer  is  given  by 


H-  1 

y{t)  =  Xm{t  -  tm),  (2.38) 

fflsO 

where  Om  represents  the  weighting  coefficients  and  tm  the  time  delays.  In  this  example 
Xm{t)  represents  the  carrier  modulated  signal  (not  the  baseband  representation  used 
throughout  the  rest  of  this  report)  which  is  given  by, 

3^(0  =  (2.39) 

where  Cmit)  represents  the  complex  amplitude  of  the  signal  measured  at  sensor  m,  tj 
represents  the  spatial  frequency,  0  represents  the  carrier  frequency,  and  Vm  represents  the 
displacement  from  sensor  0  in  a  direction  parallel  to  the  reference  baseline  (note  that 
Tm  —  md  for  a  uniform  linear  array). 

If  the  values  of  the  time  delays  tm  are  small  (i.e.  tm  «  0.51  BW  where  BW 
represents  the  bandwidth  of  the  modulating  signal  cj^t))  then  the  beamformer  output  can 
be  accurately  approximated  by, 


y{t)  =  Xm(i)  (2.40) 

m*  0 

Additionally,  by  choosing  the  time  delays  so  that, 

(2-41) 

the  beamformer  equation  becomes, 

y(0=  ^OmXJit)  e^y  (2.42) 

m=  0 


Inspection  of  this  last  result  shows  that  for  appropriate  choice  of  time  delays,  the 
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beamformer  output  y{t)  is  the  spatial  Fourier  transform  of  the  weighted  sensor  data  at  time 
instance  t  (similar  to  equation  (2.13)).  That  is,  the  beamformer  is  an  alternate 
implementation  of  the  classical  methods  described  previously. 

Following  the  same  procedure  as  before  to  compute  the  DF  spectrum  (see  section 
2.2),  then 


T”!  iV*!  iV“l 

S(4>)  =  X  [( X“"  «■’"”)] .  (2  «) 

<s0  m-Q  m=0 


or  in  the  more  general  case, 


h<i>)  -  7^  2 

j  —  rt 

L  a  -A  ^ 

(2.44) 

^  “  0 

Expressed  in  matrix  form 

m  -  0  TB  -  U 

k<t')  =  jygWg, 

(2.45) 

where  the  vector  g  is  defined  as 

g  = 


(2.46) 


and  the  matrix  Y  was  defined  in  section  2.4. 


Equation  (2.44)  is  identical  in  form  to  the  classical  spectral  estimator  described 
by  equation  (2.21)  except  that  the  windowed  data  matrix  Y  is  used  in  place  of  the  matrix 
X,  and  the  vector  g  is  a  more  general  version  of  the  steering  vector  e.  The  two  vectors 
are  identical  if  equation  (2.41)  is  satisfied  or  the  elements  of  the  vector  g  are  chosen  so 
that 


9n  — 


(2.47) 


where  P„  and  B  were  described  in  section  2.5.  In  terms  of  the  time  delay  coefficients 
used  in  the  beamforming  network,  this  is  equivalent  to. 


4.= 


(2.48) 


where  c  is  the  propagation  speed  of  light. 
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3.0  THE  ESTIMATED  AUTOCORRELATION  MATRIX 

The  autocorrelation  matrix,  which  was  first  introduced  in  conjunction  with 
classical  direction  finding  algorithms,  is  a  central  feature  in  the  derivation  of  the 
superresolution  algorithms  to  follow.  All  these  algorithms  are  derived  assuming  the  true 
autocorrelation  is  known,  which  is  rarely  true  in  practice.  The  autocorrelation  matrix  is 
normally  estimated  from  the  available  input  data,  and  the  manner  in  which  this  is  done 
can  impact  on  the  performance  of  the  superresolution  algorithm  especially  in  the  case  of  a 
small  number  of  sensors.  Given  the  importance  of  the  estimation  procedure,  some  of  the 
more  common  methods  are  discussed  here  along  with  some  of  their  drawbacks. 

Estimation  of  the  autocorrelation  matrix  using  only  a  single  data  sample  will  be 
discussed  first  since  it  simplifies  the  task  of  illustrating  some  of  the  concepts  involved. 
Estimation  involving  several  data  samples  measured  at  different  instances  in  time  is 
introduced  later,  followed  by  a  discussion  of  the  effects  of  noise. 

Before  proceeding  with  a  description  of  the  autocorrelation  matrix  estimation 
methods,  it  is  useful  to  examine  the  structure  of  the  true  autocorrelation  matrix  as  this  will 
provide  a  basis  for  comparing  various  estimation  methods.  The  definition  of  the 
autocorrelation  matrix  is  given  by, 


R=E{XX"},  (3.1) 

where  X  is  a  data  matrix  (or  vector)  and  is  described  in  the  following  sections.  The 
dimensions  of  R  are  defined  here  to  be  (p+l)*(p+l)  where  the  value  of  p  is  less  than  the 
number  of  sensors  N.  The  definition  of  p  has  been  chosen  this  way  to  correspond  to  the 
filter  order  of  the  all  pole  filter  methods  discussed  later  on  in  this  report.  If  we  assume  that 
the  received  signal  is  made  up  of  M  signals  with  unique  bearings  plus  uncorrelated  sensor 
noise  (either  spatial  and/or  temporal),  then  the  data  matrix  has  the  form 

u 

X  =  ^S,  +  N.  (3.2) 

t  si 

where  Sk  is  the  data  matrix  formed  from  the  signal  in  the  absence  of  all  other  signals  or 
noise.  Noting  that  the  noise  is  uncorrelated  and  that  signals  arriving  from  different 
directions  are  also  uncorrelated  (with  respect  to  position),  then  the  autocorrelation  matrix 
can  be  redefined  as 


M 

R  =  E{  ^S,S7}  +  E{NN“}.  (3.3) 

I  =i 

Since  the  signal  is  assumed  to  be  deterministic  (not  random)  the  signal  correlation  matrix 
is  defined  in  this  report  as, 


w 

R,  =  ^S,S«, 

I  si 


(3.4) 


and  the  noise  correlation  matrix  as 
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(3.5) 


Rf,  =  E{NN“}. 


An  important  aspect  of  the  true  autocorrelation  matrix  R  is  the  rank  of  the  signal 
correlation  matrix  R,.  If  the  number  of  signals  is  less  than  the  number  of  rows  or  columns 
in  R  (i.e.  M  <  p)  then  R,  will  have  rank  M.  Based  on  this,  it  would  be  possible  to  set  up  M 
equations  from  the  rows  or  columns  of  Rj  to  solve  for  the  bearings  of  the  M  unknown 
signals  exactly.  This  is  the  basis  for  improved  performance  of  superresolution  spectral 
estimators  compared  to  classical  estimators  which  are  based  on  correlating  an  ideal  signal 
with  the  data  which  requires  substantially  more  data. 

In  practical  situations,  the  true  autocorrelation  matrix  R  is  not  usually  available 
and  must  be  estimated  from  the  data.  In  this  case  the  bearings  cannot  be  solved  exactly, 
but  must  also  be  estimated  generally  using  least  mean  square  techniques.  The  accuracy  of 
the  bearings  then  becomes  a  function  of  the  accuracy  of  the  autocorrelation  matrix 
estimate. 

3.1  ESTIMATING  THE  AUTOCORRELATION  MATRIX 

The  estimated  autocorrelation  matrix  is  formed  from  the  input  data  using  the 
following  expression: 


R  =  XX".  (3.6) 

The  various  methods  of  estimating  the  autocorrelation  matrix  discussed  in  the  following 
sections  differ  only  in  the  manner  of  setting  up  the  data  matrix  X. 

In  the  simplest  case,  the  data  matrix  is  a  vector  defined  as 


*/  = 


xo 

ii 

X2 


(3.7) 


Xn-i 


In  this  case  the  lower  case  form  is  used  to  denote  a  vector  and  the  subscript  /has  been 
added  to  denote  the  forward  data  case  (i.e.,  X  =  nj).  In  the  forward  case,  the  data  is 
ordered  from  0  to  N-l.  In  a  similar  manner,  a  backward  data  vector  may  also  be  defined 
as, 


Xt  = 


Xn-Z 


* 

lo 


(3.8) 


The  backward  formulation  follows  from  the  form  of  the  backwards  signal  model.  For 
example,  letting  le  forward  data  vector  elements,  i/n,  be  represented  by  equation  (13), 
the  backward  elements  can  be  represented  as, 
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(3.9) 


Ijn  —  ^  Cm(  t)G  ^  Cni(  t)  C 

m-  1  Tn=  1 

In  comparing  the  forward  and  backward  signal  models,  the  only  difference  in  the  two 
expressions  are  the  signal  phases  represented  by  0m  and  Qm  respectively.  Since  the  apparent 
signal  bearings  are  unchanged,  the  bearings  determined  using  either  data  vector  will  be 
identical.  The  above  analysis  does  not  follow  for  noise,  so  that  in  the  noisy  case  the 
bearings  estimated  using  either  data  vector  will  only  be  approximately  the  same. 

One  difficulty  in  using  either  Xjt  or  x*  to  estimate  the  autocorrelation  matrix  using 
equation  (3.6)  is  that  since  only  one  unique  vector  is  used  to  make  the  estimate,  the 
resulting  matrix  will  have  rank  1.  Superresolution  algorithms  depend  on  the  rank  of  the 
signal  correlation  matrix  equaling  the  number  of  signals,  so  if  the  number  of  signals  is 
greater  than  1,  resolution  is  severely  degraded. 

A  second  difficulty  with  computing  the  estimated  autocor'elation  matrix  in  this 
manner  is  that  the  elements  of  the  matrix  are  sensitive  to  instantaneous  noise 
perturbations  in  the  data,  whereas  the  true  autocorrelation  matrix  is  only  sensitive  to  the 
the  statistics  of  the  noise  which  are  assumed  to  be  constant  or  very  slowly  changing  over 
the  measurement  period. 

Methods  to  deal  with  these  difficulties  are  discussed  in  the  following  sections. 

3.2  THE  AUTOCORRELATION  METHOD 

One  method,  called  the  autocorrelation  method  [3-1],  is  to  append  p  zeros  to  the 
start  and  p  zeros  to  the  end  of  the  data.  Variants  include  appending  zeros  to  the  start  or 
end  only.  The  data  is  then  constructed  from  a  number  of  shifted  subarrays  of  size  p-1-1 
where  p  <  Vas  shown  below. 


■  0, 

0, 

0,  . 

■•)  0>  Xq, 

...,  Xfi-  p-l,  Xn-p, 

. . . ,  Xfi-  3 , 

Xn-2,  ly-1 

0, 

0, 

0,  .. 

. . ,  Xq,  Xi, 

Xff-p,  Xf/-p*i, 

...,  Xf,.2, 

1 ) 

0 

0, 

0, 

0,  ., 

Xi,  X2, 

Xn-ptl,  Xn-pt2, 

...,  Xn-U 

0, 

0 

,(3.10) 

0, 

0, 

Xq  , 

■  Ip- 3)  Ip- 2) 

...,  Xft-Z,  Xfi-7, 

...,  0, 

0, 

n 

0, 

Xq  , 

Xi,  ., 

•I  Xp.  2,  Ip-1> 

...,  Xs-2,  Xn-1, 

...,  0, 

0, 

0 

.  lo, 

I2, 

• )  Xp.  1 ,  Ip, 

...,  Xfi- 1 ,  0 , 

...,  0, 

0, 

0 

where  the  subscript  f  is  again  used  to  denote 

the  forward  case. 

The  corresponding 

conjugate 

backward  matrix  is  given  by 
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X*_  1 


0, 

0, 

0, 

0, 

Xn-U 

. . . ,  Ip, 

Ip-  1 , 

...,  I2, 

Xl, 

xo 

0, 

0, 

0, 

...,  lAf-l, 

Xh-2, 

...,  ip-i, 

Ip- 2 1 

....  Ii, 

lo, 

0 

0, 

0, 

0, 

...,  Xn-Z, 

Xh-3, 

...,  Ip-2) 

Ip-  3, 

...,  IQ) 

0, 

0 

0. 

0, 

iv-ii 

•  ••I  Xff-  pt2t 

Xh-  p*i, 

...,  I2, 

Xl, 

...,  0, 

0, 

0 

0, 

Xh-I, 

Xn-2, 

•  •  •  1  Xn-  p*  1  > 

Xh-p, 

...,  Ii, 

lO, 

...,  0, 

0, 

0 

Xn-1, 

Xh-2, 

Xh-Z, 

...,  Itf-pj 

Xh-  p-ii 

...,  Ifl, 

0, 

...,  0, 

0, 

0 

(3.11) 


Although  the  autocorrelation  method  avoids  the  rank  deficiency  problem 
discussed  in  the  previous  section,  resolution  of  spectral  estimates  (induding  dassical 
methods)  are  degraded  due  to  the  pre-  and  postwindowing  of  the  aata  (i.e.  the  values  of  the 
unknown  data  in  the  data  matrix  n,  x.2,  x-z,  x.p  and  xh,  xntj^.i  are  assumed 

to  be  zero).  This  problem  becomes  more  severe  as  the  number  of  sensors  is  decreased. 
Figure  3.1  illustrates  one  example  of  this  where  three  equal  power  signals  with  bearings  of 
40,  50,  and  120  degrees  are  intercepted  by  an  8  dement  sensor  array  with  one  half 
wavdength  spadng.  The  signal  to  noise  ratio  was  65  dB.  Three  methods  (the 
autocorrdation  method  and  two  others  to  be  described  in  the  following  sections)  were  used 
to  estimate  the  autocorrdation  matrix  for  p  =  4.  The  Thermal  Noise  estimator  (section 
5.2.6)  was  then  used  to  compute  the  corresponding  DF  spectrums.  In  the  case  of  the 
autocorrelation  method,  the  signals  at  40  and  50  degrees  were  not  resolved  despite  the  high 
signal  to  noise  ratio.  In  general,  the  results  were  poor  compared  with  the  other 
autocorrelation  matrix  estimation  methods. 


Based  on  the  poor  resolution  of  this  method,  the  autocorrdation  method  is 
considered  too  inaccurate  for  systems  with  small  arrays  (e.g.  tactical  systems)  and  is  not 
discussed  in  the  rest  of  this  report. 


FIGURE  3.1:  DF  Spectrum  based  on  three  different  methods  of  autocorrelation 
matrix  estimation 
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3.3 


THE  COVARIANCE  METHOD 


Another  method  of  autocorrelation  matrix  estimation  is  called  the  covariance 
method  [3-1].  It  operates  only  on  the  known  data  and  thus  avoids  the  windowing  problem. 
Dividing  the  data  into  overlapping  subarrays  of  size  p+1,  (a  technique  which  is  called 
spatial  smoothing  [3-2])  the  forward  data  matrix  is  given  by 


X 


1 

f  ^N-p 


lo. 

®1> 

X2,  — , 

Xu-p-1 

Xu 

Xz, 

Xn-p 

12, 

ZZy 

I4,  •••, 

Xh-p*1 

(3.12) 

Xp, 

Xp*2t  ■•■1 

Xh-1 

The  corresponding  backward  matrix 


X*-  1 


is  given  by, 


Xs-U 

Xu-2, 

Xg-z, 

Xp 

Xn-2, 

X/f-3, 

Xg-i, 

...,  Ip-l 

Xn-Z, 

Xu- 4, 

Xg-s, 

...,  1^2 

Xn-p-U 

Xg-p-Z, 

Xg-p-Zt 

...,  Zq 

(3.13) 


In  the  covariance  method,  assuming  the  data  is  noisy  and  only  a  single  sensor 
sample  is  available,  the  rank  of  the  estimated  autocorrelation  matrix  will  be  the  smaller 
value  of  either  N-p,  which  represents  the  number  of  columns  of  the  data  matrix  X/  or  X*,  or 
p+l,  which  represents  the  number  of  rows  of  the  data  matrix.  Ideally  the  rank  will  be 
greater  than  or  equal  to  M  for  optimum  performance  of  the  superresolution  estimators. 
Conversely,  the  maximum  number  of  signal  bearings  that  may  be  estimated  is  p  subject  to 
the  constraint  that  this  value  is  less  than  or  equal  to  the  rank  of  the  data  matrix.  TUs 
constraint  can  be  expressed  as. 


p  <  N-p  and  p  <  p-Hl.  (314) 

Since  the  rightmost  expression  is  always  true,  the  largest  value  for  p  is  found  by  solving  the 
leftmost  expression.  The  corresponding  maximum  number  of  signals  (where  M  <  p)  then  is 
given  by, 


M<^.  (3.15) 

It  should  be  noted  that  spatial  smoothing  effectively  decreases  the  sensor  aperture 
from  to  p+1.  This  results  in  a  decrease  in  resolution  due  to  the  smaller  effective 
aperture,  although  this  is  partially  offset  by  the  averaging  effect  of  the  extra  columns  in  the 
data  matrix. 

An  example  of  the  performance  of  the  covariance  method  is  shown  in  Figure  3.1. 

A  complete  description  of  this  example  is  given  at  the  end  of  section  3.3. 
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3.4 


THE  MODIFIED  COVARIANCE  METHOD 


A  more  recently  developed  technique,  called  the  modified  covariance  method 
[3-3], [3-4],  doubles  the  number  of  data  vectors  used  to  form  the  data  matrix  X,  by 
combining  both  forward  and  backward  data  matrices  of  the  covariance  method.  The  result 
is  given  by. 


X/>  =  ilX;.X*l-  (316) 

where  the  subscripts  fb  are  used  to  denote  forward-backward,  and  the  data  matrices  Xf  and 
Xft  are  defined  by  equations  (3.12)  and  (3.131  respectively.  The  modified  covariance 
method  is  also  sometimes  called  the  forwara-bac^ard  method. 

In  the  modified  covariance  method,  assuming  the  data  is  noisy  and  only  a  single 
sensor  sample  is  available,  the  rank  of  the  estimated  autocorrelation  matrix  will  be  the 
smallest  value  of  either  2»{N-p),  which  represents  the  number  of  columns  of  Xfb  (double 
that  of  X/  or  Xj),  or  p+1  which  represents  the  number  of  rows  of  Xff  The  maximum 
number  of  sign^  bearings  that  may  be  estimated  is  p  subject  to  the  constraint  that  the 
rank  is  equal  to  or  greater  than  the  number  of  signals.  The  constraints  can  be  expressed  as, 

p  <  2(N-p)  and  p  <  p-l-1.  (3-17) 

Since  the  rightmost  expression  is  always  true,  the  largest  value  for  p  is  found  by  solving  the 
leftmost  expression.  The  corresponding  maximum  number  of  signals  (where  M  <  p)  then  is 
given  by, 


M<^.  (3.18) 

In  comparing  this  expression  to  the  corresponding  expression  for  the  covariance  method 
(equation  (3.15))  the  advantage  of  the  modified  covariance  method  is  clear.  A  greater 
number  of  bearings  can  be  estimated  from  a  single  sensor  sample  without  sacrificing  as 
much  resolution  due  to  the  decreased  aperture  size. 

An  example  of  the  performance  of  the  modified  covariance  method  compared  to 
the  preAriously  described  methods  is  given  in  Figure  3.1.  A  complete  description  of  this 
example  is  given  at  the  end  of  section  3.3. 

3.5  TIME  AVERAGING 

In  the  case  where  a  number  of  time  samples  are  available,  a  better  estimate  of  the 
autocorrelation  matrix  may  be  achieved  simply  by  time  averaging  as  shown  here: 

“  =  T  Z  (3-36) 

t  =0 

where  R(t)  is  the  autocorrelation  estimate  formed  for  time  sample  t.  Equivalently,  the 
data  vector  can  be  modified  in  the  following  manner. 
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(3.20) 


X  =  ^[X(0),  X(1),  X(2), X(r-I)  ], 

where  X{t)  is  the  data  matrix  (as  described  in  the  previous  sections)  formed  from  time 
sample  t.  In  this  form  the  relationship  between  time  averaging  is  obvious,  i.e.  time 
averaging  is  averaging  performed  over  time,  and  spatial  smoothing  is  averaging  performed 
over  position. 

The  advantage  of  time  averaging  is  that  the  resultant  estimated  autocorrelation 
matrix  becomes  less  sensitive  to  the  effects  of  temporal  noise.  Additionally,  in  the  case  of 
signals  that  are  uncorrelated  in  time  (i.e.  they  are  independent  signals  transmitted  from 
separate  transmitters),  averaging  increases  the  number  of  linearlv  independent  columns  of 
the  matrix  X  by  order  T  compared  to  the  submatrices  X(0j,  X(l),  etc.  In  this  case  the  rank 
deficiency  problem  can  be  overcome  simply  by  taking  a  sufficient  number  of  time  samples 
(i.e.  T  >  Musing  the  covariance  method  and  T>  Mjh  using  the  modified  covariance 
method  without  spatial  smoothing  or  p  =  N-1)  to  achieve  the  required  rank  of  the 
estimated  signal  correlation  matrix,  without  having  to  resort  to  spatial  smoothing. 

Finally,  in  the  case  of  correlated  signals  (e.g.  multipath)  the  signals  do  not 
decorrelate  in  time  so  that  spatial  smoothing  technique  must  be  used. 

3.6  THE  EFFECT  OF  NOISE  ON  AUTOCORRELATION  MATRIX 

ESTIMATION 

One  way  to  observe  the  effects  of  noise  on  the  estimation  of  the  autocorrelation 
matrix  is  to  determine  the  mean  and  variance  of  the  matrix  elements  when  they  are 
estimated  using  noisy  data.  Starting  with  the  covariance  method  with  time  averaging,  but 
no  spatial  smoothing,  the  elements  of  the  estimated  autocorrelation  matrix  can  be  defined 
as, 

* 

=  T  S  (3-21) 

<  =0 

If  the  data  is  assumed  to  be  corrupted  by  white  Gaussian  noise,  each  sensor  data  value  can 
be  represented  as  the  sum  of  a  signal  plus  a  noise  component,  that  is, 

X„{t)  =  S„{t)  +  Vn{t).  (3.22) 

Substituting  this  relationship  back  into  equation  (3.21)  gives, 

T  -  1 

r,;  =  Y  X  (3.23) 

i  =0 

The  mean  value  of  the  elements  when  *  ^  j  is  given  by, 

1  * 

E{''„}  =  (3.24) 

t  =0 


and  when  i  =  j  (the  main  diagonal  elements), 
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(3.25) 


T-l  ^ 

I  =0 


where  represents  the  noise  power. 

The  elemental  variance  is  given  by, 

u  =  E{  |ri^-E{r4y}l^}  =  (3-26) 

/  *0  J  *0 

The  above  result  is  based  on  the  fact  that  the  variance  of  a  process  VY*  is  o-y*  and  the 
variance  of  the  process  Jf  y  is  (Tx^(Ty^  where  X  and  Y  represent  uncorrelated  white  Gaussian 
processes  with  variances  and  respectively.  To  simplify  the  above  variance  expression 
a  new  parameter,  k,j,  is  defined  so  that, 

T-  I 

~  5^  )>  (3.27) 

t  *0 

where  s^  is  the  sum  of  the  individual  signal  powers.  (Note  that  for  a  large  number  of 
samples,  and  uncorrelated  signals,  Kij »  1).  Using  this  definition  then,  the  variance 
becomes. 


.  =  .  (3.28) 

Since  bearing  accuracy  is  a  function  of  the  ratio  of  the  input  noise  power  (a^)  to 
signal  power  (s^),  equation  (3.28)  can  be  normalized  by  dividing  through  by  s^  (which  is  the 
equivalent  maximum  "sign^  power"  of  the  autocorrelation  elements)  and  reexpressing  the 
normalized  variance  in  terms  of  the  input  signal  to  noise  ratio.  Calling  this  the  normc^zed 
elemental  variance  for  the  covariance  meth(M,  the  result  is  given  by 

Vc  =  ](2KijSNRr^  +  SNR:^,  (3.29) 

where  the  input  signal  to  noise  ratio  measured  at  a  sensor  is  given  by 

SNR  =  1^.  (3.30) 

Note  that  Vi^  is  a  measure  of  the  signal  to  noise  power  ratio  of  the  elements  of  the 
autocorrelation  matrix  (as  opposed  to  the  signal  to  noise  power  ratio  of  the  data). 

Inspection  of  equation  (3.29)  shows  that  the  normalized  variance  is  an  inverse 
function  of  SNR  for  signal  to  noise  ratios  greater  than  zero  and  an  inverse  function  of  SNR^ 
for  signal  to  noise  ratios  less  than  zero.  Figure  3.2  illustrates  this  effect  through  simulation 
of  a  5  element  array  with  half  wavelength  spacing,  p  =  4,  and  T  =  5.  The  bearing  of  the 
incoming  signals  was  40,  50,  and  120  degrees,  and  they  were  uncorrelated.  The  variance 
was  calculated  from  1000  simulation  runs  performed  for  input  signal  to  noise  ratios  ranging 
from  -60  to  -1-60  dB  in  1  dB  steps  and  averaged  for  all  the  elements  of  the  autocorrelation 
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matrix. 


The  significance  of  this  effect  is  that  the  bearing  accuracy  will  degrade  more 
rapidly  for  signd  to  noise  ratios  less  than  zero,  indicating  the  need  to  remove  as  much  noise 
as  possible  from  the  data  before  estimating  the  autocorrdation  matrix. 


FIGURE  3.2:  Elemental  variance  as  a  function  of  signal  to  noise  ratio 


In  the  more  general  case  where  spatial  smoothing  and/or  time  averaging  is  used 
equation  (3.29)  can  easily  be  modified  to  become, 

Vc  =  ^2KijSNIt^  +  SNR:^).  (3.31) 

Here  K  represents  the  total  number  of  terms  averaged  together  and  can  be  expressed 
mathematically  as 


K=T\N-v).  (3.32) 

If  the  modified  covariance  method  is  used  twice  as  many  terms  are  involved  in  the 
estimation  of  the  autocorrelation  matrix.  For  comparison  purposes  it  is  useful  to  keep  the 
same  expression  for  K  and  modify  equation  (3.31)  instead.  At  first  sight  this  suggests  a 
simple  relationship  between  the  elemental  variance  for  the  covariance  method  (v^  and  the 
elemental  variance  for  the  modified  covariance  method  (Vm),  namely, 


=  O.SVc  (3.33) 

This  expression  is  not  always  valid  as  explained  in  the  following  analysis. 

A  closer  inspection  of  the  modified  covariance  method  reveals  that  it  can  be 
defined  in  terms  of  the  covariance  method  ^^.s, 
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(3.34) 


—  0.5( Cij  +  Cmn)) 

where  m  =  p-j  and  n  =  p-i  In  this  last  expression  rjj  represents  an  element  of  the 
estimated  autocorrelation  matrix  determined  using  the  modified  covariance  method  and  c,j 
represents  an  element  determined  using  the  covariance  method.  Since  each  element  of  the 
modified  covariance  estimate  is  the  average  of  two  elements  of  the  covariance  estimate,  the 
modified  covariance  estimate  would  be  expected  to  have  half  the  variance  (as  predicted  by 
equation  (3.33))  as  long  as  the  errors  in  Cij  and  Cmn  are  uncorrelated. 

There  are  two  conditions  where  the  elements  are  correlated.  The  first  case  is  for 
elements  lying  on  the  main  cross  diagonal  of  the  autocorrelation  matrix  (i.e.  rop,  rM, 
r5p-2,  Tpo).  In  this  case  Cij  and  Cmn  represent  the  same  elements  so  that  equation  (3.34) 
simplifies  to 


r,j  —  Ctj.  (3.35) 

Consequently  for  the  cross  diagonal  elements,  the  normalized  variance  is  actually  given  by 
equation  (3.31). 

The  second  case  where  the  elements  are  correlated  occurs  when  spatial  smoothing 
is  used  in  conjunction  with  the  modified  covariance  method.  Due  to  the  overlapping  nature 
of  the  subarrays,  the  elements  Cij  and  Cmn  are  formed  from  many  of  the  same  terms  (where 
each  term  has  the  form  i^x/).  By  comparing  how  the  two  elements  Cj,  and  c^n  are  formed, 
an  expression  for  the  number  of  common  terms  can  be  derived,  namely,  hT  where 

(N-p-\p-i-j\  if  the  result  is  >  0 
h  =  (3.36) 

I  0  otherwise 


Since  the  averaging  operation  in  equation  (3.34)  has  no  effect  on  the  hT  common  terms, 
then  the  variance  will  be  Vc  while  the  improvement  due  to  the  remaining  K-  hT 
uncorrelated  terms  will  be  0.5Vc.  From  this,  the  elemental  variance  for  the  modified 
covariance  method  can  be  defined  as. 


_  VchT  +  0.5  ViiK-  hT) 
K 


or  in  terms  of  SNR, 


II.  =  ^(1  +  57^)  {2K,fiNR'  +  SNR\ 


(3.37) 


(3.38) 


If  the  value  of  /i  =  0,  then  equation  (3.33)  applies.  If  on  the  other  hand  h  >  0  then 
equation  (3.38)  can  be  rewritten  as. 


Vm  =  2^ (2  -  {2Ki^NR^  +  SNR:^  for  h  >  0.  (3.39) 

Inspection  of  this  result  shows  that  for  large  values  of  N  [i.e.  N  »  2p),  the  normalized 
elemental  variance  of  the  modified  covariance  method  degrades  to  that  of  the  covariance 
method. 
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Given  that  the  elemental  variances  are  not  all  necessarily  equal,  it  is  useful  to 
determine  the  average  elemental  variance  since  this  quantity  provides  a  more  useful 
measure  of  how  well  the  estimated  autocorrelation  matrix  approximates  the  true 
autocorrelation  matrix.  To  do  this,  the  assumption  is  made  that  =  1.  For  a  large 
number  of  trials  involving  signals  with  unifor^y  distributed  phases  this  assumption  is 
reasonable.  For  a  single  trial  involving  a  very  limited  number  of  sensor  samples  and/or 
correlated  signals,  the  formulas  derived  in  the  following  discussion  will  only  be 
approximations. 

In  the  case  of  the  covariance  method  the  elemental  variances  are  all  equal  so  that, 

T,  =  (3.40) 

where  the  overbar  is  used  to  denote  the  mean  value.  The  situation  is  not  as  straight 
forward  for  the  modified  covariance  method  since  the  value  of  h  in  equation  (3.38)  changes 
with  each  element.  After  some  algebraic  manipulations,  however,  the  result  is  given  by, 

-  ^  (  4  +  3JV  +  forp>f  {3,41) 

and 

Vm  =  (1  -  6(A'-^p)t/+  ij  ^OT  P<J  (3.42) 

Again,  as  N  increases  for  a  fixed  subarray  size,  or  fixed  value  of  p,  the  elemental  variance  of 
the  modified  covariance  estimate  approaches  that  of  the  covariance  method. 


FIGURE  3.3:  Blow  up  of  the  elemental  variance  shown  in  Figure  3.2 
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For  the  example  shown  in  Figure  3.2,  the  predicted  value  of  Vm  —  0.6Uc  based  on 
equation  (3.41)  which  is  in  excellent  agreement  with  the  simulated  results  (see  Figure  3.3). 
For  signal  to  noise  ratios  much  greater  than  0,  this  translates  into  an  equivalent  increase  in 
the  signal  to  noise  ratio  of  2.2  dB  when  using  the  modified  covariance  method  compared  to 
the  covariance  method. 

The  concepts  embodied  by  equations  (3.40),  (3.41),  and  (3.42)  are  illustrated  in 
Figures  3.4  and  3.5.  These  figures  show  the  decrease  in  the  variance  as  a  function  of  the 
number  of  terms  for  time  averaging  (Figure  3.4)  and  spatial  smoothing  (Figure  3.5)  when 
either  the  covariance  or  modified  covariance  methods  were  used.  In  the  time  averaging 
case,  the  data  samples  were  assumed  to  be  taken  from  a  5  element  array  with  half 
wavelength  spacing.  Each  sensor  sample  was  assumed  to  be  uncorrelated  with  the  previous 
sample.  In  the  spatial  smoothing  case  the  samples  were  taken  as  overlapping  subarrays  (5 
elements)  from  a  single  snapshot  of  a  very  large  array.  The  bearings  of  the  incoming  signals 
were  40,  50,  and  120  degrees  and  they  were  all  uncorrelated.  Statistics  were  computed  from 
1000  trials  for  each  value  of  K. 

As  predicted  by  equation  (3.40),  the  decrease  in  variance  for  the  covariance 
method  is  a  function  of  the  factor  IjK  whether  time  averaging  or  spatial  smoothing  was 
used.  The  same  result  holds  true  for  the  modified  covariance  method  when  time  averaging 
is  used  as  predicted  by  equations  (3.41)  and  (3.42).  In  the  spatial  smoothing  case  the 
advantage  of  the  modified  covariance  method  over  the  covariance  method  tegins  to 
disappear  as  K  increases,  exactly  as  predicted  by  these  equations.  The  theoretical  results 
for  Figures  3.4  and  3.5  based  on  the  theoretical  equations  are  not  shown  since  they  were 
virtually  indistinguishable  from  the  simulated  results. 
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- 


Number  of  Terms  Averaged  (K) 


FIGURE  3.5:  Normalized  variance  using  spatial  smoothing 


The  results  presented  here  do  not  predict  the  ultimate  bearing  accuracies  of  any 
particular  superresolution  DF  method,  since  the  bearing  estimation  procedure  is  highly 
nonlinear  (although  linear  approximations  are  possible  at  high  signal  to  noise  ratios). 
However,  these  results  are  useful  for  predicting  some  of  the  ways  in  which  noise  affects 
bearing  estimation  as  well  as  providing  a  measure  of  the  merits  of  various  autocorrelation 
matrix  estimation  methods. 
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4.0  THE  MODELLING  APPROACH 


The  key  to  improving  the  performance  of  DF  estimators  over  that  of  classical 
methods  lies  in  taking  better  advantage  of  the  form  of  the  sensor  data.  As  discussed  in 
section  2,  determining  the  bearing  of  received  signals  is  equivalent  to  the  problem  of 
determining  the  frequencies  of  complex  sinusoids  in  noise.  Since  the  characteristics  of  the 
signal  and  noise  are  different,  they  can  be  modelled  separately. 

4.1  THE  SIGNAL  MODEL 

For  a  single  signal  in  a  noiseless  environment,  the  data  ffom  the  sensor  in  an 
N  sensor  system  can  be  represented  by  equation  (1.1)  which  is  repeated  here  as 


Xn=  (4.1) 

An  alternate  representation  for  this  equation  is  given  by. 


Zn  —  (XoXn-i, 


(4.2) 


where  the  coefficient 


(4.3) 


and  is  easily  computed  from  the  data.  The  advantage  of  this  alternate  representation  of 
the  data  is  that,  at  least  in  this  case,  it  provides  a  simple  method  of  extending  the  data 
sequence  and  corresponding  autocorrelation  sequence  indefinitely.  The  Fourier  transform  of 
the  infinitely  extended  autocorrelation  sequence  then  results  in  the  ideal  DF  spectrum. 
Although  the  single  bearing  could  also  be  determined  directly  from  oo,  the  problem 
becomes  more  difficult  when  several  signals  are  involved. 

In  the  multiple  signal  environment,  again  assuming  no  noise,  the  sensor  data  can 
be  represented  by. 


x,»=  ,  (4.4) 

m»  I 

where  the  subscript  m  is  used  to  distinguish  between  the  M  signals.  The  equivalent 
alternate  representation  in  this  case  is  given  by, 

M 

3^  —  ^^Un*Xn.m-  (^•®) 

ms  1 

The  relationship  between  the  coefficients  represented  by  Om  and  the  spatial  frequencies  of 
the  signals  is  not  as  clear  cut  as  in  the  single  signal  case.  However,  a  relationship  does  exist 
as  demonstrated  in  the  following  analysis. 

To  simplify  this  analysis,  equation  (4.4)  is  rewritten  as. 
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where  for  simplicity 


M 

^  ~  ^^CmPm  > 
m*  1 


(4.6) 


Cm  —  Cnj(t), 

and  the  complex  signal  poles  pm  are  defined  as 

*}Umd 

Pm=  e 


(4.7) 


(4.8) 


Since  the  complex  amplitude  c„  is  of  no  interest  for  bearing  determination,  and  noting 
equation  (4.6)  forms  a  linear  set  of  equations  given  by, 


Xo  - 

Cl 

C2 

C3 

-  CAf  = 

0 

-1 

-1 

-1 

-  1 

0 

Xi  - 

CiPi  - 

C2P2  - 

C3P3  -  ••• 

-  CuPm  = 

Xj  - 

-2 

ClPl  - 

-2 

C2P2  - 

-2 

C3P3  -  ... 

-  2 

-  CmPu  = 

0 

(4.9) 

-AT*! 

-.V*l 

-Af»l 

-y*i 

0 

Xff-l  - 

ClPl  - 

C2P2  - 

C3P3  -  ... 

-  CuPnf  = 

then  Cm  can  be  eliminated  using  standard  techniques.  For  example,  to  remove  the  first 
coefficient  Ci  from  anv  row  (where  each  of  the  above  equations  is  referred  to  as  a  row  and 
are  ordered  as  shown),  the  row  is  multiplied  by  pi  and  then  subtracted  from  the  previous 


row.  If  the  operation  is  performed  on  the  last  N-l  rows,  the  following  set  of  N-1  equations 
results: 

(xo  -  iipi) 

-  C2(l  -Pipi*) 

-C3(l-piP3‘) 

...  -cj<(l  -pxp'u)  =  0 

(ii  -  3;2Pi) 

-  C2ip2  -piP2^) 

-  C3(P3*  -piP3^) 

•••  -  c«(pm‘ -piPv*)  =  0 

{X2-X3P1) 

-  2  -  3 

-  C2(P2  -  PlP2  ) 

-  2  *3 

-  C3(P3  -piP3  ) 

-  2  -  3 

•••  -Ct^pn  -piPu  )  =  0  (4.10) 

{Xft.2  -  Xs-lPl) 

.  -y»2  -y*!. 

-  C2(P2  -  P1P2  ) 

,  -y*2  -y^u 

-  C3(P3  -  PiP3  )  - 

,  -y42  -y*!. 

•••  -  Ct^pu  -  PiPu  )  =  0 

Removing  any  of  the  other  coefficients  proceeds  in  an  identical  manner.  That  is,  to  remove 
Ck  from  a  row,  multiply  the  row  by  p*  and  subtract  it  from  the  previous  row.  Note  that 
each  time  this  operation  is  performed,  the  resultant  set  of  equations  is  reduced  by  1 
equation. 


If  the  procedure  outlined  above  is  carried  through  until  all  the  coefficients  Cm  have 
been  removed,  then  the  resultant  N-M  equations  have  the  form  given  by  equation  (4.5) 
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where 


=  Pi  +  P‘2  +  Pi  +  ■■■+  Pu 
0-2  =  -P\P2  -  PlPi  -  •••  -  P2Pi  -  —  -  PU-\PU 

03  =  P1P2P3  +  P1P2P4  +  •••  +  P2P3P4  +  ...  +  Plt.2Pu-iPu  (4.11) 


au-l  =  (-1)"(P2P3P4...PM  +  PlP3P4...Plf  +  ...  +  PiP2P4...Pjm) 
dU  =  (-1)*^*' (PiP2P3P4P5  -  Pm) 

As  in  the  single  signal  case,  once  the  coefficients  represented  by  a„  have  been  determined, 
the  data  sequence  i„  can  be  extended  indefinitely.  The  resultant  DF  spectrum  (computed 
from  the  extended  data  set  using  either  the  direct  or  indirect  methods  discussed  in  section 
1)  can  be  used  to  exactly  determine  the  signal  bearings.  Figure  4.1  shows  an  example  of  the 
improvement  in  the  DF  spectrum  using  tMs  technique  compared  to  classical  methods. 


FIGURE  4.1:  Comparison  of  the  DF  spectrum  generated  using  the  Linear 

Prediction  method  and  the  Classical  method  for  two  signals  of 
equal  power  at  40  and  50  degrees. 
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4.2  THE  NOISE  MODEL 

In  the  case  where  noise  is  present,  but  no  signals,  the  sensor  output  may  be 
represented  by 


Xn=T}r,{t).  (4.12) 

Since  noise  is  not  deterministic  (i.e.,  not  completely  predictable),  it  must  be  handled 
statistically.  For  example,  for  complex  white  Gaussian  noise,  the  autocorrelation  sequence 
is  given  by. 


and 


r„(m)  =  0  for  m  ^  0, 


(4.13) 


r„(0)  =  (4.14) 

where  is  the  variance  of  the  noise  process.  This  model  is  useful  for  modelling  internal 
sensor  noise  which  is  typically  white  Gaussian  ncise  (in  the  temporal  sense)  with  the  same 
variance  but  uncorrelated  between  sensors.  It  is  also  useful  for  modelling  external  (e.g. 
atmospheric  noise)  omnidirectional  noise  with  equal  power  in  all  directions. 

For  diffuse  external  noise  sources  which  have  an  unequal  noise  distribution  with 
direction,  the  noise  can  be  modelled  as  filtered  complex  white  Gaussian  noise.  That  is, 

K 

In  “  ^  ' bjnVn-mi 
in  =  0 

where  bo  =  1,  K  represents  the  order  of  the  noise  process,  and  Vk  represents  a  complex 
Gaussian  white  noise  process  with  a  variance  Since  Vk  is  not  a  deterministic  signal,  but  a 
stochastic  process,  some  estimation  method  must  be  used  to  determine  the  optimum  value 
of  the  coefficients  (as  opposed  to  determining  the  exact  values  of  a,n  for  the  signal  only 
case  discussed  in  the  previous  section). 

The  choice  of  the  filter  order  Kin  equation  (4.15)  generally  depends  either  on  the 
known  characteristics  of  the  noise  (e.g.  for  white  noise  K  =  0),  or  is  limited  by  the  amoimt 
of  available  data.  Since  the  autocorrelation  sequence  can  only  be  estimated  to  lag  N,  then 
K<  N-l. 

4.3  THE  SIGNAL  PLUS  NOISE  MODEL 

One  approach  to  improving  spatial  frequency,  which  follows  from  the  previous 
discussion  in  the  preceding  sections,  is  to  combine  the  signal  and  noise  models  to  give  a 
model  capab'e  of  handling  the  signals  plus  noise  problem.  One  such  model  is  called  an 
autoregressive  moving  average  (ARMA)  filter  and  is  given  by. 


In 


Ir*-m  4"  ^ 


m~  1 


m- 0 


(4.16) 


Ideally  the  autoregressive  (AR)  filter  coefficients  represented  by  Om  are  chosen  to  model  the 
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signals,  and  the  moving  average  (MA)  filter  coefficients  represented  by  bm  are  chosen  to 
model  the  noise.  In  practice,  this  is  not  always  true,  or  possible.  The  methods  used  to 
estimate  these  values  are  the  basis  of  the  various  superresolution  DF  estimators  discussed 
in  this  report. 

4.4  COMPUTING  THE  DF  SPECTRUM 

Once  the  values  for  Om  and  bm  have  been  estimated,  the  data  sequence,  and 
correspondingly  the  autocorrelation  sequence  can  be  extended  indefinitely  by  computing 
the  ur^nown  values  of  Xn.  The  Fourier  transform  of  the  extended  autocorrelation  sequence 
then  gives  the  power  spectral  density  function. 

Computationally,  the  direct  method  of  computing  the  power  spectral  density 
function  from  the  data  sequence  is  simpler  and  is  given  by 

=  X{<I>)X{<I>)*  (4.17) 

where  A((/))  is  the  Fourier  transform  of  the  extended  data  sequence. 

X{<j>)  may  also  be  computed  from  equation  (4.16)  by  taking  the  Fourier  transform 
of  both  sides  to  give, 


Xi<p)  =  -^am  X{<P)  +  ^bm  V(<P) 

m  s 1  m~0 

and  then  rearranging  to  get 

where 


A{<p)  =  1  +  '^Om  e^'^, 

m=  1 


and 


m  = 

m=  0 


(4.18) 

(4.19) 

(4.20) 

(4.21) 


Substituting  equations  (4.3)  and  (4.5)  back  into  equation  (4.2),  the  power  spectral 
density  function  can  be  computed  based  on  the  coefficients,  a,,  and  That  is. 


(4.22) 


Since  for  a  white  noise  process 


31 


(4.23) 


where  is  the  variance  of  the  noise,  then 

S(0)  =  . 

The  matrix  representation  of  equation  (4.24)  is, 


(4.24) 


S{<l>)  =  2(d 


e!?bb"e, 

e^"ep’ 


(4.25) 


where  ep  and  e,  are  p+1  and  ^+1  element  steering  vectors  (described  in  section  2.3.1 
and  defined  by  equation  (2.18)),  the  autoregressive  filter  coefficient  vector  a  is  defined  as 


n 


ai 


and  the  moving  average  coefficient  vector  b  is  defined  as 


f  6o 

bi 


(4.26) 


(4.27) 


For  spectral  estimation  purposes,  equation  (4.25)  is  useful.  However,  since 
ultimately  the  goal  is  to  determine  signal  bearings,  a  simpler  form  of  the  DF  spectrum  can 
be  used.  For  example,  the  main  interest  in  the  DF  spectrum  is  its  shape  (i.e.  to  locate  the 
simal  peaks),  and  consequently  only  the  relative  values  of  the  actual  spectrum  are  needed. 
Therefore  the  noise  coefficient  2a^  can  be  ignored  resulting  in  the  expression. 


(4.28) 


A  further  simpbfication  can  be  made  based  on  the  observation  that  if  the 
coefficients  represented  by  bm  are  chosen  to  model  the  noise  only,  they  provide  no 
information  on  the  location  of  the  signal  peaks  and  in  fact,  could  make  it  more  difficult  to 
determine  the  location  of  the  true  peaks  by  masking  them.  This  observation  does  not 
simplify  the  task  of  calculating  these  coefficients,  but  it  does  result  in  the  simplified  DF 
spectrum  given  by, 


1 

e^*^p 


(4.29) 
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5.0  DF  ESTIMATORS 


In  the  following  sections,  DF  estimators  which  are  inherently  based  on  the 
modelling  concepts  discussed  in  section  4,  are  discussed.  In  filter  terminology  these 
estimators  may  be  divided  into  three  classes,  namely,  all  zero,  all  pole,  and  pole-zero  [3-1]. 

5.1  ALL  ZERO  ESTIMATORS 

In  all  zero  estimators,  or  more  commonly  called  moving  average  fMA)  estimators, 
only  the  moving  average  part  of  the  ARM  A  filter  given  by  equation  (4.16)  model  is  used. 
The  values  of  the  autoregressive  parameters  a„  are  set  to  equal  0,  giving, 

n  =0 

where  v*  represents  a  complex  white  Gaussian  noise  source.  The  DF  spectrum  can  be 
derived  from  equation  (4.28)  by  noting  that  the  coefficients  a„  all  equal  to  zero  resulting  in 
the  expression. 


S{<(>)  =  e^bV  (5.2) 

From  the  discussion  in  section  4,  the  MA  model  was  shown  to  be  appropriate  for 
modelling  noise-like  processes  in  the  DF  spectrum  (e.g.  spatial  noise  which  has  broad 
spectral  peaks  and  sharp  nulls).  Although  it  was  also  shown  in  section  4  that  signals  can  be 
accurately  modelled  using  an  all  pole  filter  model,  the  MA  model  can  also  be  used  to  model 
signals,  albeit  with  reduced  efficiency.  That  is,  a  large  number  of  filter  coefficients, 
compared  to  the  numbers  of  signals,  may  be  required  to  provide  an  accurate  DF  spectrum. 

Further  insight  into  the  properties  of  MA  estimators  can  be  gained  comparing 
equation  (5.2)  to  the  classical  estimator  defined  by, 

^  r„(n)e'^‘^,  (5.3) 

n*  -  q 

where  estimation  of  the  autocorrelation  lags  is  discussed  in  section  1.  Provided  that  the 
autocorrelation  sequence  results  in  a  positive  spectrum  {S{<t>)  >  0  for  all  (p)  then  equation 
(5.3)  can  be  factored  into  the  form  given  by  (see  also  Appendix  A), 


S(«  =  .  (5.4) 

n  =  0  n  =  0 

where  in  this  case  the  coefficients  represented  by  b„  are  computed  from  the  autocorrelation 
lags.  The  matrix  form  of  this  expression  is  identical  to  equation  (5.2).  In  other  words, 
although  the  underlying  development  philosophy  is  different,  classical  DF  estimators  (as 
described  in  this  report)  are  a  subclass  of  moving  average  estimators  [5-1]. 

From  this  analysis,  it  is  apparent  that  the  performance  of  moving  average 
estimators  would  not  be  expected  to  significantly  improve  on  the  performance  of  classical 
estimators.  In  general,  the  performance  limitations  of  moving  average  estimators  can  be 
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viewed  as  a  failure  of  the  model  to  extend  the  data  sequence  beyond  the  known  data.  A 
consequence  of  this  fact  is  that  the  white  noise  process  v*  in  equation  (5.1)  is  not 
predictable,  and  so  unknown  values  cannot  be  predicted. 

5.2  ALL  POLE  ESTIMATORS 

In  the  all  pole  model,  only  the  autoregressive  part  of  the  ARMA  filter  defined  by 
equation  (4.16)  is  used.  The  values  of  the  moving  average  coefficients  b„  are  set  to  0, 
giving, 


n  =  1 

From  the  discussion  in  section  4,  it  is  apparent  that  this  model  is  more  appropriate  for 
generating  DF  spectra  which  contain  signal  peaks,  than  are  MA  techniques.  As  a  result,  the 
^1  pole  estimators  discussed  in  the  following  sections  are  generally  superior,  in  terms  of 
accuracy,  for  direction  finding  purposes  than  are  MA  and  classical  techniques.  As  a  result, 
these  estimators  are  often  called  superresolution  DF  methods.  The  differences  in  the 
following  methods  are  in  the  manner  that  the  filter  coefficients  On  are  selected  (although  in 
some  of  these  methods  On  is  not  calculated  directly). 

5.2.1  Autoregressive  Method 

The  Autoregressive  (AR)  method  is  based  on  defining  the  autocorrelation  sequence 
using  equation  (5.5)  and  the  relationships  defined  in  section  2.1.1  to  give, 


3,(m)  =  J’rj(^-Ji)  for  all  m>0, 

n  =  I 


(5.6) 


and  for  m=  0, 


(5-7) 

n  ■  I 

Equations  (5.6)  and  (5.7)  are  known  as  the  Yule-Walker  equations  or  normal  equations, 
and  are  also  sometimes  referred  to  as  the  discrete-time  Wiener-Hopf  equations.  Once  the 
values  of  the  coefficients  a„  have  been  determined,  the  autocorrelation  sequence  can  be 
extended  infinitely,  and  an  improved  estimate  of  the  spectrum  calculated. 

Equations  (5.6)  and  (5.7)  can  be  incorporated  into  a  single  matrix  equation,  called 
the  augmented  normal  equations,  to  give, 

Ra  =  a\,  (5.8) 

where  R  is  the  (p-f-l)*(p-l-l)  augmented  autocorrelation  matrix,  the  coefficient  vector  a 
was  defined  previously  oy  equation  (4.26),  and  u  is  a  p+1  element  unit  vector  defined  as. 
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u  = 


(5.9) 


1 
0 
0 

6 

Assuming  the  autocorrelation  matrix  R  is  invertible  (in  the  presence  of  white  Gaussian 
noise  R  will  be  full  rank  and  invertible  although  when  estimated  the  result  may  not  be), 
equation  (5.8)  can  be  rewritten  in  terms  of  a  as, 


a  =  o^R-y  (5.10) 

In  cases  where  R  is  not  invertible,  the  method  of  solution  of  a  is  discussed  in  section  7.1. 

Alternatively,  if  only  the  location  of  spectral  peaks  and  their  power  with  respect  to 
the  rest  of  the  spectrum  is  required,  it  is  only  necessary  to  solve  for  the  coefficients  a„ 
using  the  linear  set  of  equations  represented  by  equation  (5.6).  In  matrix  form  this  set  of 
equations  can  be  expressed  by. 


RpW  =  -r, 


(5.11) 


where  Rp  is  the  p»p  normal  autocorrelation  matrix,  w  is  the  coefficient  vector  defined  by 


and  r  is  the  vector  defined  by 


tti 

02 


(5.12) 


(5.13) 


Assuming  the  autocorrelation  matrix  Rp  is  invertible,  then  equation  (5.11)  can  be  rewritten 
in  terms  of  w  as 


w  =  (5.14) 

In  cases  where  Rp  is  not  invertible,  the  method  of  solving  w  is  described  in  section  7.2. 

The  vector  r  can  also  be  defined  by  noting  that  the  augmented  autocorrelation 
matrix  R  can  be  partitioned  in  terms  of  r  and  Rp  as 
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R  = 


rii(0) :  rH 


:  R 


■p  J 


(5.15) 


Although  these  two  definitions  of  r  are  equivalent,  the  second  definition  is  more  useful 
when  the  estimated  autocorrelation  matrix  (discussed  in  section  3)  is  used  in  place  of  the 
true  autocorrelation  matrix.  In  this  case  the  true  values  in  equation  (5.15)  are  simply 
replaced  by  their  appropriate  estimates. 


A  simple  relationship  also  exists  between  the  coefficient  vectors  a  and  w,  namely, 


a  = 


(5.16) 


Once  the  coefficients  have  been  determined,  the  power  spectral  density  function 
can  be  calculated  using  equation  (4.29)  which  is  repeated  here  as. 


In  terms  of  the  augmented  autocorrelation  matrix  and  the  solution  for  the  vector  a  given  in 
equation  (5.9),  the  DF  spectrum  for  the  Autoregressive  method  may  also  be  expressed  as. 


Sar{<P)  = 


(5.18) 


where  the  scaling  factor  has  been  ignored. 

5.2.2  Maximum  Entropy  Method 

The  Maximum  Entropy  (ME)  method  [5-2]  is  closely  related  to  the  Autoregressive 
method.  In  this  method  the  extrapolation  of  the  autocorrelation  sequence  is  made  in  such  a 
way  as  to  maximize  the  entropy  of  the  data  series  represented  by  the  sequence.  The  data 
series  would  then  be  the  most  random,  in  an  entropy  sense,  of  all  possible  series  which 
include  the  known  autocorrelation  lags  as  part  of  the  sequence. 

The  entropy  rate  for  a  Gaussian  random  process  is  proportional  to 


ln[  S{<p)  ]  du,  (5.19) 

where  the  spatial  frequency  u  is  defined  in  terms  of  the  bearing  (p  by  equation  (1.2), 

Wo  =  it/ d,  and  the  power  spectral  density  function  is  represented  by 


S(0)  =  ^  rj*(m) 


(5.20) 


To  maximize  the  entropy,  the  derivative  of  equation  (5.19)  is  taken  with  respect  to  the 
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unknown  autocorrelation  lags  (i.e.  r„{m)  where  \m\  >  p  are  the  unknown 
autocorrelation  lags).  This  leads  to 


du  =  0  for  I  m|  >  p. 


(5.21) 


Equation  (5.21)  implies  that 


has  a  finite  Fourier  expansion,  that  is, 


(5.22) 


where  Cm  =  c^m-  The  summation  term  on  the  right  hand  side  of  this  expression  can  be 
factored  (see  Appendix  A),  as  long  as  5(0)  >  0  for  all  0,  and  the  resultant  expression 
inverted  to  give. 


m  = 


(5.23) 


The  matrix  form  is  given  by, 


Sui{4>) 


(5.24) 


From  this  analysis,  it  is  clear  that  the  maximum  entropy  method  belongs  to  the 
class  of  all  pole  DF  estimators.  Additionally  it  has  been  shown  [5-3]  that  for  the  problem  of 
signals  in  white  Gaussian  noise,  and  a  uniform  linear  antenna  array,  the  Maximum 
Entropy  method  is  identical  to  the  Autoregressive  method  (i.e.  5jir^0)  =  S/uj(0)).  For 
other  types  of  noise  or  antenna  spacings,  the  two  methods  are  not  identical. 

In  the  case  of  non-uniform  antenna  spacing,  the  maximum  entropy  solution  for  the 
set  of  equations  represented  by  equation  (5.21)  usually  requires  some  form  of  gradient 
search  technique.  This  can  lead  to  a  number  of  practical  difficulties  which  are  not 
addressed  in  this  report. 


5.2.3  LLaear  Prediction  Method 


In  time  series  modelling  the  Linear  Prediction  (LP)  method  predicts  either  a 
future  or  past  data  value  using  a  sequence  of  current  data  values.  In  Erection  finding,  this 
is  equivalent  to  predicting  either  the  first  sensor  (backward  prediction),  the  last  sensor 
(forward  prediction),  or  both  the  first  and  last  sensor  (forward-backward  prediction)  in  a 
group  of  sensors.  These  three  types  of  predictors  are  described  in  the  sections  5.2.3. 1  to 
5.2.3.3. 
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5.2.3. 1  Forward  Linear  Prediction 


Mathematically,  the  forward  prediction  case  can  be  expressed  as, 


fi»l 


(5.25) 


where  the  subscript  /  is  used  to  denote  the  forward  prediction  case.  The  error  in  the 
estimate  of  Xm  is  given  by. 


efm=Xm-Xm.  (5.26) 

The  values  of  the  coefficients,  af„,  ate  determined  by  minimizing  the  variance  of 
the  error  given  by, 

Vf=E{\ef,n\'}  =  E{ej,,e%},  (5.27) 

where  the  error  is  assumed  to  be  a  zero  mean  process.  To  minimize  this  value,  the 
derivative  is  taken  with  respect  to  each  of  the  coefficients,  ajk,  giving 


=  «■ 


(5.28) 


which,  using  the  results  from  equations  (5.25)'^5.27)  and  given  that  =  0,  simplifies  to 


o/* 


B{ef„  aw*}  =  0. 


(5.29) 


Replacing  Cfm  by  equation  (5.26),  Xm  by  equation  (5.25),  and  expanding  gives, 

E{xn  aw*}  +  E{  ^o/„  aWn  aw*}  =  0,  (5.30) 

»*  I 

where  0  <  ib  <  p.  This  results  in  a  system  of  linear  equations  that  can  be  expressed  in  terms 
of  the  autocorrelation  parameters  as, 

rii(k)  +  ^a/n  rj^k-n)  =  0.  (5.31) 

n  s  I 

Additionally,  by  incorporating  equation  (5.26)  into  the  right  side  of  equation 
(5.27)  and  expanding,  the  optimum  variance  is  given  by 


*♦ 

Replacing  Xm  by  the  conjugate  of  equation  (5.25)  and  then  applying  the  result  from 
equation  (5.29)  in  order  to  simpbfy,  then 
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(5.33) 


Again  expanding  in  terms  of  equation  (5.26), 


—  ^m}  "f"  ^  ^Q/n 


(5.34) 


and  re-expiessing  the  result  in  terms  of  the  autocorrelation  parameters, 


=  r„(0)  +^a/n  r4-n). 


(5.35) 


The  system  of  equations  described  by  equations  (5.31)  and  (5.35)  can  be 
represented  in  matrix  form  as, 


(5.36) 


which  is  identical  in  form  to  equation  (5.8),  the  Autoregressive  model  system  equations.  In 
fact,  for  the  problem  of  signals  in  white  Gaussian  noise,  the  Linear  Prediction  method  and 
the  Autoregressive  method  are  identical  [5-41.  The  corresponding  DF  spectrum  then  is 
given  by  equation  (5.18)  (i.e.  Sip(<p)  =  S^/t(^)). 

In  the  case  where  the  true  autocorrelation  matrix  is  unknown,  the  matrix  R  in 
equation  (5.36)  is  replaced  by  the  estimated  autocorrelation  defined  by. 


R  =  X,X/ 


(5.37) 


where  X/  is  the  forward  data  matrix  given  by  one  of  equations  (3.7),  (3.10),  or  (3.12)  (the 
concept  of  forward  and  backward  data  originated  in  linear  prediction  research 
[5-2], [5-5115-6]). 

5.2.3.2  Backward  Linear  Prediction 

The  case  of  the  backward  prediction  coefficients  proceeds  in  much  the  same 
manner,  where  the  backward  estimate  is  defined  as 


flftn  Xmtn, 


(5.38) 


and  the  subscript  b  is  used  to  denote  the  backward  prediction  case.  The  final  result  is 
given  by. 


x{-k)  +  ^fl6n  r„{n-k)  =  0  for  0<A<p, 


(5.39) 
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for  k=0. 


(5.40) 


atm  r*x(n)  =  Vt,^^ 


n  -  1 


The  equivalent  matrix  solution  then  is, 


is, 


T 

Ra6  = 


(5.41) 


By  taking  advantage  of  the  fact  that  the  autocorrelation  matrix  is  Hermitian,  that 


r’’  =  R*,  (5.42) 

and  complex  conjugating  both  sides  of  equation  (5.41)  gives 


Since  this  is  identical  in  form  to  the  Autoregressive  equations  (5.8),  the  expression  for  DF 
spectrum  for  the  backward  linear  predictor  is  the  same  as  for  tne  Autoregressive  method 
given  by  equation  (5.18). 

Equation  (5.43)  is  also  identical  in  form  to  equation  (5.36),  so  that  the  solution  for 
the  backward  coefficients  and  error  variance  can  be  expressed  in  terms  of  the  forward 
values  as. 


and 

an4=o*/.  (5.45) 

As  in  the  forward  case,  when  the  true  autocorrelation  matrix  is  unknown,  an 
estimate  is  used  in  its  place  which  is  defined  by. 


R  =  XiXj,  (5.46) 

where  Xt  is  the  backward  data  matrix  given  by  one  of  equations  (3.8),  (3.11),  or  (3.13). 
5.2.3.3  Forward-^Backwaid  Linear  Prediction 

In  practice,  when  the  estimated  autocorrelation  matrix  is  used,  the  conjugate 
relationship  between  forward  and  backward  coefficients  expressed  by  equation  (5.45)  is  not 
true  if  the  coefficients  are  calculated  separately.  However,  by  constraining  these 
coefficients  to  satisfy  this  relationship. 


Onj  —  flni  —  flni 

and  solving  the  forward  and  backward  prediction  equations  (5.J6)  and  (5.43) 
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simultaneously,  a  better  estimate  of  the  coefficients,  for  spectral  estimation  purposes, 
results.  This  method  is  called  the  Forward-Backward  Linear  Prediction  (FBLPj  method. 

The  solution  can  be  derived  by  minimizing  the  quantity  Vj  +  Vi,  subject  to  the 
constraint  given  by  equation  (5.47).  Following  the  same  sort  of  derivation  as  used 
previously  in  the  forward  or  backward  cases  the  result  in  matrix  form  is, 

Ra  =  (5.48) 

Again  the  DF  spectrum  is  given  by  equation  (5.18)  (i.e.  Sii{<p)  =  SAn{(f>)). 

The  estimated  autocorrelation  matrix  used  in  the  least  mean  square  solution  is  the 
modified  covariance  estimate  defined  by, 


R  =  X/Xj  +  XftX!  =  X/jX/i, 


(5.49) 


where  X^j  is  given  by  equation  (3.16). 

5.2.4  Minimum  Variance  Method 

The  Minimum  Variance  (MV)  method  is  based  on  the  output  of  a  beamformer 
which  passes  all  energy  arriving  from  the  look  direction  and  adaptively  minimizes,  in  an 
optimal  manner,  the  energy  arriving  from  all  other  directions  [5-7].  This  is  equivalent  to 
minimizing  the  variance  of  the  beamformer  output  subject  to  the  constraint  that 
beamformer  gain  is  unity  in  the  look  direction. 

The  mathematical  derivation  proceeds  as  follows.  Consider  the  output  of  a 
beamformer  given  by. 


J/n 


**0 


(5.50) 


where  there  are  N  antennas  in  the  antenna  array,  and  each  output  y„  is  formed  from  a 
subarray  of  p  antennas.  In  matrix  notation  the  system  of  equations  embodied  by  equation 
(5.50)  can  be  rewritten  as, 

y  =  Xc,  (5.51) 

where  the  output  vector  y  is  defined  by 


y  = 


7^ 


yo 

y\ 

y2 

yn-p 


(5.52) 


the  data  matrix  X  is  identical  to  the  forward  data  matrix  X/  defined  by  equation  (3.12) 
and  the  beamformer  weights  vector  c  is  given  by. 
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The  beamformer  output  variance  is  given  by 

V  = 

which  can  be  rewritten  in  terms  of  the  input  data  as, 


(5.53) 


(5.54) 


V  =  E{c“XX"c}  =  c“Rc,  (5.55) 

where  R  is  the  autocorrelation  matrix. 

Assuming  the  antenna  array  is  uniform  and  linear,  the  transfer  function  of  the 
beamformer  is  given  by, 


H{(f>)  =  ^  c*  (5.56) 

*s0 

Given  that  the  direction  of  interest  is  represented  by  (pa,  then  the  beamformer  gain  will  be 
constrained  so  that 


H{(po)  =  1. 


(5.67) 


Expressed  in  matrix  form. 


e“(^o)c  =  c"e(^o)  =  1,  (5.58) 

where  e(^)  is  the  p+1  element  steering  vector  e  (defined  by  equation  2.19)  evaluated  at 
<(>  =  <p0‘ 


The  beamformer  coefficients,  Ck,  are  then  derived  by  minimizing  the  variance 
defined  in  equation  (5.55)  subject  to  the  constraint  represented  by  equation  (5.58).  The 
solution  technique  for  tWs  problem  is  to  use  the  Lagrange  multiplier  [5-8]  which  involves 
minimizing  the  expression 


F  =  e*Rc  +  S{e*e(<po)  -  1),  (5.59) 

with  respect  to  the  coefficient  c*’  (minimizing  with  respect  to  S  yields  the  original 
constraint  equation  given  by  equation  (5.58)).  Performing  this  minimization  yields  the 
system  of  equations  expressed  in  matrix  form  as, 

Rc  +  5e((^o)  =  Q,  (5.60) 

where  0  is  a  p+1  element  column  vector  whose  elements  are  all  0. 


42 


The  final  solution  for  the  beamformei  coefficients  comes  by  first  solving  equations 
(5.58)  and  (5.60)  to  eliminate  c  and  then  expressing  6  in  terms  of  e(0o)  and  R,  which  , 
gives, 

6  = - — - ,  (5.61) 

e”(0o)R‘'e(^o) 

then  substituting  this  expression  back  into  equation  (5.60)  and  solving  for  c  gives 


e(^o)V‘e(4^o)’ 


(5.62) 


Finally,  the  minimum  variance  equation  can  be  rewritten  by  incorporating  this 
expression  for  c  back  into  equation  (5.55).  The  result  is. 


(5.C3) 


This  equation  represents  the  power  output  of  the  beamformer  for  the  chosen  look  direction 
corresponding  to  the  spatial  frequency  4>o-  The  beamformer  power  output  for  any 
direction,  then,  is  given  by, 


e  R  e 

An  important  difference  between  this  estimator  and  the  other  estimators  discussed 
is  that  the  Minimum  Variance  estimator  determines  the  power  (using  equation  5.64)  of  the 
signal  in  the  look  direction,  not  the  power  density.  The  advantage  is  that  the  power  of  the 
signals  can  be  determined  by  the  height  of  the  peaks  in  the  DF  spectrum.  The  disadvantage 
is  poorer  resolution  (see  section  5.4). 

The  all  pole  nature  of  this  estimator  is  easily  shown  by  noting  that  the 
denominator  term  in  equation  (5.64)  can  be  factored  as  (see  Appendix  A), 

e“R'*e  =  e“aa"e  (5.65) 


In  the  case  where  an  estimate  of  R  is  used,  a  necessary  condition  is  e“R'‘e  >  0  for  all  <f> 
(see  Appendix  A).  The  resultant  estimator  then  has  the  all  pole  form  represented  by 
equation  (4.29). 

5.2.5  Thermal  Noise  Method 

The  Thermal  Noise  method  (TN)  is  based  on  a  beamformer  which  functions  in  a 
similar  way  to  the  Minimum  Variance  beamformer  except  the  gain  is  not  constrained  to 
unity  in  the  look  direction  [5-9].  For  example,  based  on  the  previous  derivation  of  the 
Minimum  Variance  methoa,  the  weight  co^dent  vector  can  be  defined  in  terms  of 
equation  (5.60)  to  get. 


c  =  -fiR  ‘e, 


(5.66) 
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where  6  was  originally  defined  by  equation  (5.61).  If  the  gain  constraint  for  the  look 
direction  is  chosen  to  be  some  arbitrary  nonzero  function  with  respect  to  the  look  angle  0o, 
equation  (5.66)  remains  unchanged,  but  6  will  be  modified  according  to  the  new  gain 
constraint.  In  other  words,  the  purpose  of  6  is  to  control  the  gain  of  the  array  in  the  look 
direction. 

In  the  Thermal  Noise  method  the  parameter  6  is  set  to  -1  with  the  resulting  weight 
coefficient  vector  given  by. 


c  =  Re.  (5.67) 

For  this  choice  of  weight  vector  the  beamformer  adaptively  minimizes  all  energy  arriving 
from  all  directions.  The  output  in  this  case  is  noise  only,  that  is,  in  equation  (5.50)  the 
beamformer  output  represents  a  noise  process  which  has  been  called  "thermal  noise". 

Since  this  approach  is  very  similar  in  concept  to  the  Autoregressive  Method,  the 
DF  spectrum  is  computed  in  an  almost  identical  manner.  That  is. 


S(«  =  — 1 — . 


(5.68) 


where  in  this  case  is  the  overall  beamformer  transfer  function  (compared  to  equation 
(5.56)  which  is  the  transfer  function  for  a  particular  look  angle).  The  function  H[(f)yH{(f)) 
can  be  determined  by  deriving  the  output  power  spectrum  in  response  to  a  white  noise 
input  with  a  variance  of  1.  For  example,  given  the  input  white  noise  process  U{<p)  the 
output  spectrum  is  given  by. 


V{<P)  =  (5.69) 

and  the  output  power  spectrum  is  given  by, 

(5.70) 

* 

Since  t/(0)  is  a  white  noise  process  with  variance  of  1,  then  U{<l>)  U{<p)  =  1.  Simplifying 
equation  (5.70)  then, 


V{4>)*  V[<t>)  =  (5.71) 

In  other  words  for  a  white  noise  input,  the  output  spectrum  is  identical  to  the  beamformer 
transfer  function. 

The  matrix  form  output  power  for  equation  (5.71)  is  given  by  equatior  (5.55) 
which  is  repeated  here  as, 


V(0)*V(^)  =  c”R,c.  (5.72) 

Using  the  fact  that  the  autocorrelation  matrix  for  a  white  noise  process  is  given  by  a^l 
where  the  variance  <72=1^  then  equation  (5.72)  simplifies  to. 
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V(^)*V(^)  =  c"c.  (5.73) 

Finally  using  the  results  from  equations  (5.67),  (5.68),  (5.71),  and  (5.73),  the  DF  spectrum 
for  the  Thermal  Noise  method  is  given  by, 

=  (5.74) 

eK  e 

As  in  the  case  of  the  MV  estimator,  the  all  pole  nature  of  this  estimator  is  easily 
shown  by  noting  that  the  denominator  term  in  equation  (5.74)  can  be  factored  as. 


e“R‘^e  =  e"aa'fe,  (5.75) 

which  assumes  that  S((^)  >  0  for  all  </>  (see  Appendix  A).  The  resultant  estimator  then  has 
the  all  pole  form  represented  by  equation  (4.29). 


5.3  POLE-ZERO  ESTIMATORS 

In  section  4.3  it  was  shown  that  the  autoregressive  moving  average  (ARMA)  filter 
was  best  suited  for  modelling  signal  in  noise  problems.  Ideally  the  moving  average  (or  all 
zero)  part  of  the  filter  models  the  noise  and  the  autoregressive  (or  all  pole)  part  models  the 
signal.  In  practice  development  of  computationally  fast  algorithms  based  on  the  ideal 
ARMA  have  not  been  as  successful  as  the  all  pole  models.  Current  methods  typically  rely 
on  using  time  consuming  search  algorithms,  which  are  largely  inappropriate  for  real  time 
applications  (at  least  until  faster  hardware  is  available).  For  this  reason,  these  types  of 
estimators  are  not  considered  in  this  report. 

One  pole-zero  based  method,  the  Adaptive  Angular  Response  method,  is 
considered  since,  as  will  be  seen,  this  method  is  a  simple  modification  of  an  all  pole 
method. 

5.3.1  Adaptive  Angular  Response 

The  Minimum  Variance  method  can  be  modified  to  give  the  power  density 
(instead  of  power)  by  dividing  the  Minimum  Variance  beamformer  output  by  the  effective 
noise  beamwidth  of  the  beamformer  to  get  the  average  power  density  in  the  beam.  In  the 
Adaptive  Angular  Response  (AAR)  method,  this  technique  is  used  to  give  an  estimate  of 
the  true  spectral  power  density  function  [5-10], [5-11].  Expressed  mathematically, 

(5.76) 

where  Suv{(p)  was  defined  previously  by  equation  (5.64),  and  Bjv(0)  is  the  effective  noise 
bandwidth  of  the  beamformer. 

Evaluated  for  a  particular  look  direction  (i.e.  0  =  0o)  the  effective  beamwidth  can 
be  calculated  using. 
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f  \n(.4‘)\' i'P 

0a)  =  — - 5 — 


(5.77) 


where  the  bearing  angle  <p  is  expressed  in  radians.  The  transfer  function,  E{<p),  for  a 
Minimum  Variance  beamformer  was  previously  defined  by  equation  (5.56),  and  can  be 
expressed  in  matrix  form  as, 

H{<f>)^  ce.  (5.78) 

Incorporating  this  expression  into  equation  (5.77)  and  recalling  that  H{(f>o)  =  1,  then 


H  H  H 

<Po)  =  e  CoCoe  d<f>  =  coCo, 
Jo 


(5.79) 


where  co  is  the  coefficient  vector  c  evaluated  at  (^  =  (j>Q.  Generalizing  this  3xpression  for 
any  value  of  0  and  then  substituting  back  into  equation  (5.76)  gives. 


m 


(5.80) 


Using  the  expression  for  S»v{(p)  given  by  equation  (5.64)  and  the  generalized  form  (i.e. 
any  look  direction)  of  equation  (5.62)  for  c,  the  final  result  is, 


Saar{<P)  = 


e“R'‘e 


(5.81) 


e”R-2e 

5.4  A  COMPARISON  OF  DF  ESTIMATORS 

In  this  section,  five  different  DF  methods  are  compared  which  are  representative  of 
the  approaches  discussed  up  to  this  point.  The  DF  estimators  associated  with  these 
methods  are  summarized  here  as: 


Bartlett  (section  2.3.2): 

SBAja{<P)  =  6*^116 

(5.82) 

Autoregressive  (section  5.2.1): 

Maximum  Entropy  (section  5.2.2): 

Linear  Prediction  (section  5.2.3): 

Sli{<P')  —  u  .j  u  -1 

eR.  unit  e 

(5.83) 

Minimum  Variance  (section  5.2.4): 

5w((/>)  =  u  .1 

eR.  e 

(5.84) 
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Thermal  Noise  (section  5.2.5): 


1 


(5.85) 


3tn{<I>)  = 


U  * 

5.  Adaptive  Angular  Response  (section  5.3.2):  Saar{<P)  =  (5.86) 

e  K*  e 

Note  that  for  p  <  N-1,  the  Bartlett  method  becomes  the  Welch  method  [5-12]. 

In  the  simulation  examples  that  follow  (unless  otherwise  indicated),  the  signal 
environment  consisted  of  three  signals  of  equal  j^wers  and  bearings  of  40,  50,  and  120 
degrees.  The  direction  finding  array  consisted  of  8  colinear  sensors  with  half  wavelength 
spacing.  Noise  between  sensors  was  uncorrelated.  Signal  phases  and  sensor  noise  were  also 
uncorrelated  from  trial  to  trial  where  each  trial  consisted  of  estimating  the  bearings  as  the 
signal  to  noise  ratios  were  varied  Rom  5  to  65  dB  in  1  dB  steps  (the  noise  was  scaled 
accordingly  for  each  step).  Bearing  error  variance  statistics  were  based  on  300  such  trials. 

Bearing  errors  were  calculated  by  determining  the  bearing  of  the  three  largest 
peaks  in  the  spectrum  and  subtracting  these  values  from  the  corresponding  true  values. 
Bearing  accuracy  was  then  calculated  as  the  root  mean  squared  value  of  these  errors. 
Justification  for  choosing  the  three  largest  peaks  in  the  spectrum  is  based  on  the  fact  that 
the  heights  of  the  peaks  are  generally  related  to  the  square  of  the  corresponding  signal 
power  (except  for  the  Bartlett  and  Minimum  Variance  methods  where  the  peaks  are 
proportional)  [5-4], [5-9].  Although  this  is  not  necessarily  the  optimal  method  for  choosing 
the  correct  signal  peaks,  it  is  useful  for  comparing  various  methods  and  highlighting  the 
problems  involved. 

Since  multipath  is  considered  to  be  the  most  significant  problem  in  tactical  DF  at 
VHF/UHF,  the  autocorrelation  matrix  estimates  were  generated  from  single  sample 
samples  (i.e.  T  =  1)  using  spatial  smoothing.  The  choice  of  the  model  order  p  for  spatial 
smoothing  purposes  was  mainly  based  on  choosing  the  model  order  for  which  estimator 
accuracy  was  the  highest  for  a  given  signal  to  noise  ratio  as  shown  in  Figure  5.1.  The 
variance  of  the  bearing  estimates  for  each  value  of  p  were  computed  from  300  trials  at  a 
fixed  signal  to  noise  ratio  of  50  dB. 


Model  Order  (p) 

FIGURE  5.1:  Bearing  error  variance  as  a  function  of  the  model  order  p 
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For  the  Linear  Prediction,  and  Minimum  Variance  methods  the  optimum  choice 
based  on  Figure  5.1  was  p  =  4,  and  for  the  Adaptive  Angular  Response  method  the 
optimum  choice  was  p  =  3.  Although  the  optimum  choice  for  the  Thermal  Noise  method 
using  Figure  5.1  is  p  =  4,  it  was  found  in  simulations  that  this  method  performed  better  at 
lower  signal  to  noise  ratios  for  a  model  order  of  p  =  3  with  only  a  marginal  decrease  in 
accuracy.  For  this  reason,  a  model  order  of  p  =  3  was  used  in  simulations  involving  the 
Thermal  Noise  method.  (Note  model  order  selection  for  situations  where  the  true  signal 
beairings  are  unknown  is  discussed  in  section  10). 


25r 


Direction  (degrees) 

FIGURE  5.2:  DF  Spectrum  of  the  Bartlett  method  for  the  noiseless  case 
(dashed  lines  show  the  true  signal  bearings). 


Figure  5.2  provides  am  example  of  the  DF  spectrum  generated  using  the  Bairtlett 
estimator  in  a  noiseless  environment.  As  mentioned  before,  the  resolution  of  classical  and 
moving  average  methods  is  poor  as  illustrated  in  this  example  (i.e.  the  bearings  at  40  and 
50  degrees  aire  unresolved).  For  this  reason,  the  Bartlett  estimator  is  not  considered  in  the 
following  simulation  examples. 

Of  the  matrix  estimation  schemes,  the  results  in  section  3.6  implied  that  under  the 
same  conditions,  the  modified  covariance  method  was  superior  to  the  covariance  method. 
These  results  are  confirmed  in  Figure  5.3  which  illustrates  the  improvement  in  the  Thermal 
Noise  estimator  accuracy  versus  signal  to  noise  ratio  using  the  two  different  methods.  In 
this  particular  example,  the  c'  variance  method  produced  results  which  were  about  10  dB 
poorer  in  terms  of  the  signal  to  noise  ratio  than  the  modified  covariance  method.  This  is 
worse  than  would  be  predicted  from  the  results  given  in  section  3.6  and  probably  due  to  the 
fact  that  the  autocorrelation  matrix  estimate  for  the  covariance  method  was  not  full  rank 
which  degrades  the  results  somewhat.  Results  using  other  estimators  are  similar.  For  these 
reasons,  all  further  comparisons  between  estimators  in  this  and  later  sections  is  done  using 
the  modified  covariance  method. 
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FIGURE  5.3:  Comparison  of  the  Covariance  and  Modified  Covariance  methods 
used  in  conjunction  with  the  Minimum  Variance  method. 


Figure  5.4  compares  the  performance  of  4  DF  estimators.  From  the  results  shown 
it  is  clear  that  at  high  signal  to  noise  ratios  the  estimators  are  well  behaved  with  the 
bearing  error  variance  inversely  proportional  to  the  signal  to  noise  ratio.  In  terms  of  the 
mean  elemental  variance  of  the  estimated  autocorrelation  matrix  (described  in  section  3.6), 
the  bearing  variance  is  directly  proportional  to  the  elemental  variance.  This  is  in  keeping 
with  the  comment  that  at  high  signal  to  noise  ratios  estimator  performance  will  be 
approximately  linear. 


49 


In  terms  of  estimator  accuracy,  at  high  signal  to  noise  ratios,  the  Linear  Prediction 
and  Minimum  Variance  methods  had  the  best  performance.  At  lower  signal  to  noise  ratios, 
estimator  performance  in  Figure  5.4  departs  dramatically  from  their  linear  behaviour.  The 
point  at  which  this  failure  occurs  is  called  threshold.  In  terms  of  threshold  effects,  the  best 
performance  was  achieved  equally  by  the  Adaptive  Angular  Response,  Linear  Prediction, 
and  Thermal  Noise  methods,  while  the  Minimum  Variance  method  had  the  poorest 
performance.  In  all  of  these  methods  for  this  simulation,  the  threshold  effect  was  caused  by 
merging  of  the  signal  peaks  at  40  and  50  degrees  in  the  DF  spectrum.  Consequently  the 
lower  the  threshold,  the  better  the  resolution  of  the  method. 

Although  the  Adaptive  Angular  Response  method  performed  reasonably  well  in 
these  simulations,  this  method  was  extremely  sensitive  to  model  order  as  Figure  5.1 
indicates.  At  higher  model  orders  (p  >  3)  the  performance  of  this  method  was  extremely 
poor  (e.g.  see  Figure  5.1).  This  is  due  to  two  effects,  spurious  peaks  in  the  DF  spectrum 
which  are  mistaken  for  true  signal  peaks,  and  spectral  peak  inversion.  The  spectral  peak 
inversion  problem  is  illustrated  in  Figure  5.5.  To  understand  the  cause,  it  is  useful  to  note 
that  the  spectrum  of  the  Adaptive  Angular  Response  estimator  is  simply  the  ratio  of  the 
Thermal  Noise  DF  spectrum  divided  by  the  Minimum  Variance  DF  spectrum  at  any 
particular  bearing.  Normally  the  peaks  in  the  Thermal  Noise  spectrum  are  larger  than  the 
corresponding  peaks  in  the  Minimum  Variance  spectrum  so  that  the  resultant  Adaptive 
Angular  Response  spectrum  has  positive  peaks.  Occasionally  the  opposite  is  true  with  the 
result  that  the  resultant  peaks  are  inverted,  i.e.,  valleys  are  formed  where  the  peaks  should 
be.  Given  these  problems,  the  choice  of  model  order  for  the  Adaptive  Angular  Response 
method  appears  to  be  limited  to  p  =  M  which  ensures  that  a  maximum  of  only  M  peaks 
will  exist  in  the  DF  spectrum  (no  spurious  peaks),  and  at  least  for  the  simulations 
performed  for  Figure  5.4,  avoids  the  spectrad  inversion  problem. 


FIGURE  5.5:  Spectral  peak  inversion  in  the  Adaptive  Angular  Response  DF 
spectrum  (dashed  lines  show  the  true  signal  bearings) 
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The  form  of  the  Thermal  Noise  estimator  (equation  (5.85))  provides  a  useful  basis 
of  comparison  with  the  other  estimators.  For  example  the  only  difference  between  the 
Thermal  Noise  estimator  and  the  Minimum  Variance  estimator  (equation  5.84)  is  the 
power  to  which  the  inverse  autocorrelation  matrix  is  raised.  Squaring  has  the  effect  of 
making  the  peaks  in  the  DF  spectrum  more  pronounced  which  potentially  improves  the 
resolution  capabilities  of  the  method.  This  explains  the  poor  performance  of  the  Minimum 
Variance  method  for  p  =  3  in  Figure  5.1  where  in  a  number  of  trials  the  signals  at  40  and 
50  degrees  were  not  resolved.  Raising  the  inverse  autocorrelation  matrix  to  higher  powers  is 
possible,  but  this  also  has  the  effect  of  emphasizing  spurious  peaks  which  can  degrade 
performance  of  the  estimator.  Figure  5.6  illustrates  these  effects  for  various  powers  using 
the  estimator  defined  by 


S(«  = 


-m 

e  R  e 


(5.87) 


where  m  =  1  for  the  Minimum  Variance  method  and  m  =  2  for  the  Thermal  Noise  method. 
The  generated  spectrums  have  also  been  offset  for  clarity. 


FIGURE  5.6:  Spectrum  of  the  DF  estimator  defined  by  equation  (5.87)  using 
various  power  of  m  (dashed  lines  show  the  true  signal  bearing^. 


As  mentioned  earlier,  the  choice  of  model  order  for  the  Thermal  Noise  method  in 
Figure  5.4  was  p  =  3,  even  though  accuracy  was  better  for  p  =  4  based  on  the  results 
shown  in  Figure  5.1.  It  was  found,  however,  that  in  simulations  where  p  =  4  the  threshold 
of  the  Thermal  Noise  method  degraded  to  that  of  the  Minimum  Variance  method  shown  in 
Figure  5.4.  At  higher  orders  the  accuracy  decreased  significantly  and  spurious  estimates 
were  a  problem. 

In  comparing  the  Linear  Prediction  method  to  the  Thermal  Noise  method,  it  is 
useful  to  rewrite  equation  (5.85)  for  the  Thermal  Noise  estimator  as, 
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(5.88) 


Stn{<I>)  — 


e  (  >  R  n,aiR  )e 

i  =0 


where  u*  is  a  p+1  element  column  vector  of  zeros  with  a  1  in  the  position.  In  this  form 
the  Thermal  Noise  method  is  clearly  related  to  the  Linear  Prediction  method  (equation 
(5.83)1  except  whereas  the  Linear  Prediction  method  is  based  on  using  a  single  prediction 
filter  (i  =  0  for  forward  prediction  and  t  =  p  for  backward  prediction)  the  Thermal  Noise 
method  combines  the  power  outputs  firom  p+1  filters  where  the  filter  predicts  (or 
interpolates)  the  element  in  the  sensor  subarray  being  processed.  The  outputs  &om  each 
of  these  filters  is  also  inversely  weighted  according  to  the  prediction  error  variance,  i.e.,  the 
better  the  filter  model  fit  with  the  data  the  greater  the  weighting. 

In  equation  (5.88)  each  of  the  p  filters  can  be  determined  independently.  Equation 
(5.85)  then  performs  an  average  of  the  weighted  filter  outputs.  This  is  different  from  the 
idea  of  forward-backward  linear  prediction  which  constrains  the  forward  and  backward 
filter  coefficients  to  be  related  before  the  coefficients  are  calculated. 


For  comparison  purposes  it  is  also  useful  to  define  the  null  spectrum  as, 

0(4>)  =  —  ,  (5.89) 

m 

so-called  since  nulls  in  D((p)  correspond  to  peaks  in  the  DF  spectrum.  The  null  spectrum  of 
the  Thermal  Noise  method  can  also  be  interpreted  as  the  average  of  the  null  spectrums  of 
the  p+l  individual  prediction/interpolation  filters.  A  disadvantage  of  the  Thermal  Noise 
method  is  that  in  cases  where  the  estimated  autocorrelation  matrix  is  used,  each 
prediction/interpolation  filter  will  produce  nulls  in  its  own  null  spectrum  at  slightly 
different  bearings.  Averaging  the  nulls  together,  unless  they  are  exactly  aligned,  decreases 
the  overall  null  depth  resulting  in  smaller  peaks  in  the  DF  spectrum  and  poorer  resolution 
compared  to  the  case  where  o^y  a  single  filter  is  used  (i.e.  linear  prediction).  An 
advantage,  however,  is  that  averaging  also  increases  the  immunity  to  spurious  peaks  since 
the  corresponding  spurious  nulls  in  the  individual  null  spectrums  often  don’t  occur  at  the 
same  bearing  and  as  a  result  the  corresponding  spumous  peaks  are  smoothed  out  in  the  DF 
spectrum. 

Ideally  averaging  the  null  spectrums  in  the  Thermal  Noise  method  should  also  lead 
to  more  stable  estimates.  In  practice  it  has  been  found  that  the  error  variance  and  spurious 
performance  of  the  DF  estimator  defined  by. 


5(0)  = 


1 

Hn-l  V*  ’ 
e  R  u^UiR  e 


(5.90) 


is  a  nonlinear  function  of  1 1  -  0.5p| .  That  is,  the  accuracy  and  suppression  of  spurious  peaks 
is  best  for  the  linear  prediction  case  where  the  outside  data  values  in  the  sensor  subarray 
are  being  predicted  (i  =  0  or  i=  p)  and  is  worst  when  the  middle  data  value(s)  of  the 
sensor  subarray  are  being  interpolated  {i  =  0.5p  or  i  =  0.5(p+l)).  As  a  result,  the 
accuracy  of  the  Thermal  Noise  method  was  slightly  poorer  than  the  Linear  Prediction 
method. 
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An  example  of  these  effects  is  shown  in  Figure  5.7  where  the  DF  spectrum  of  the 
Thermal  Noise  method  is  compared  to  the  spectrums  generated  using  equation  5.90.  Each 
of  the  DF  spectrums  have  been  offset  for  clarity,  and  the  spectrums  corresponding  to  z  =  3, 
and  i  =  4  have  not  been  plotted  since  these  are  identical  to  the  spectrums  corresponding  to 
i  =  0  and  i  =  1  respectively.  In  comparing  the  various  spectrums,  the  spectrum  for  z  =  2 
has  the  poorest  accuracy  (i.e.  the  peaks  corresponding  to  the  40  and  50  degree  signal 
bearings  exhibit  greater  error)  and  also  contains  a  large  spurious  peak  at  110  degrees 
(although  this  has  little  effect  on  the  Thermal  Noise  spectrum).  In  comparison,  the 
spectrum  for  i  =  0,  the  Linear  Prediction  spectrum,  has  the  best  accuracy. 


FIGURE  5.7;  DF  spectrums  for  the  Thermal  Noise  method  and  the 

corresponding  spectrums  of  the  prediction/interpolation  filters 
(dashed  lines  show  the  true  signal  bearings). 


Based  on  the  results  of  comparisons  between  the  various  DF  estimators,  two 
criteria  are  clearly  important  in  determining  estimator  performance.  The  first  is  estimator 
accuracy  above  threshold,  and  the  second  is  the  signal  to  noise  ratio  where  threshold 
occurs.  In  terms  of  accuracy,  methods  which  operate  using  higher  model  orders  achieve 
better  accuracy  above  threshold.  The  signal  to  noise  ratio  at  which  threshold  occurs  is  a 
function  of  estimator  resolution  -  the  better  the  resolving  abilities  of  the  estimator  the 
lower  the  threshold  (where  two  signals  close  in  bearing  are  considered  resolved  when  two 
corresponding  closely  spaced  peaks  are  formed  in  the  DF  spectrum).  In  terms  of  this 
criteria,  the  Linear  Prediction  method  had  the  best  performance  for  the  simulations 
summarized  by  Figure  5.4. 
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6.0  SIGNAL  AND  NOISE  SUBSPACE  DIVISION 


In  many  practical  situations,  the  sensor  data  is  the  measurement  of  a  signal  or 
signals  corrupted  by  noise.  Since  the  signal  spectrum  is  of  primary  interest,  the  addition  of 
noise  only  clouds  the  issue.  The  all  pole  type  spectral  estimators  discussed  to  this  point 
(which  include  AR,  LP,  ME,  MV,  and  TN)  are  based  on  models  which  do  not  specifically 
address  this  issue,  so  consequently,  as  the  input  signal  to  noise  ratio  decreases,  the  spectral 
estimator  performance  degrades. 

An  example  of  this  can  be  illustrated  by  considering  the  problem  of  direction 
finding  on  several  point  sources  when  the  sensor  data  is  corrupted  by  uncorrelated  complex 
white  Gaussian  noise.  The  resultant  sensor  data  can  be  represented  by. 


Xm  —  Sfn  "I"  Tlnj.  (^•^) 

Here  represents  the  uncorrupted  sensor  inputs  (i.e.  the  sensor  inputs  in  the  absence  of 
noise)  which  can  be  modelled  exactly  by  an  all  pole  model  (see  Section  4.1)  given  as. 


»  -  1 


(6.2) 


The  process  Um  represents  uncorrelated  complex  white  Gaussian  noise  with  variance 
The  exact  spectrum  of  Xm  is  given  by, 


5(*)  = - +  * 

C(^)C*(^) 

or  alternatively, 

5(0)  _  1  +  'n^C{<p)C*{<p) 

where, 

C{<p)  =  I  + 


(6.3) 


(6.4) 


(6.5) 


Clearly,  from  the  form  of  equation  (6.4),  S(<^)  is  the  spectrum  of  an  ARMA  process.  That, 
is,  the  addition  of  complex  white  Gaussian  noise  to  an  all  pole  process  results  in  an  ARMA 
process. 


From  the  preceding  analysis  it  is  apparent  that  the  ARMA  spectral  estimator 
would  be  the  most  appropriate  choice  for  estimating  the  spectrum  of  im-  However,  as 
mentioned  before,  tlus  results  in  a  difficult  non-linear  search  problem,  so  all  pole  estimators 
are  often  used  instead.  Using  any  of  the  all  pole  methods,  the  estimated  spectrum  has  the 
form, 


54 


1 


(6.6) 


k^)=  , - 

In  comparing  the  expression  for  the  estimated  spectrum  to  the  expression  for  the  true 
spectrum,  equation  (6.4),  it  is  quite  clear  that  the  estimated  spectrum  will  only  approach 
the  true  spectrum  as  the  additive  noise  power,  approaches  0.  Conversely,  the  addition 
of  noise  degrades  the  estimated  spectrum. 

If  a  method  could  be  developed  which  would  remove  most  of  the  noise  before 
processing  with  an  all  pole  spectral  estimator,  the  estimated  DF  spectrum  would  be  greatly 
improved.  A  relatively  simple  approach  is  based  on  the  observation  that  the 
autocorrelation  sequence  of  the  sensor  data  described  by  equation  (6.1),  is  given  by 

rzxim)  =  (6.7) 

Expanding  in  terms  of  Sm  and  n^, 

r„(m)  =  E{(s„.m  +  n„*m)(5„  +  n„)*}.  (6.8) 

Since  Sm  and  n,„  are  uncorrelated  then, 

r-.(m)  =  jElSn.m  3^}  +  B{n„,m  (6.9) 

Based  on  this  last  equation,  the  autocorrelation  matrix  can  be  expressed  as, 

R=R,+  R„,.  (6.10) 

where  R,  is  the  signal  autocorrelation  matrix  and  Rn  is  the  noise  autocorrelation  matrix. 
In  the  special  case  where  represents  an  additive  white  noise  process  with  a  variance 
then  the  autocorrelation  sequence  can  be  expressed  as. 


r„(0)  =  E{3„  s*}  +  772, 

(6.11) 

which  represents  the  signal  power  plus  the  noise  power,  and 

r„(m)  =  E{s„,^  s*}  for  m  #  0. 

(6.12) 

The  matrix  equation  (6.10)  becomes, 

R  =  R,  +  7/21, 

(6.13) 

where  I  is  the  identity  matrix. 

From  these  results,  the  effects  of  the  noise  can  be  removed  and  the  DF  spectrum 
estimate  enhanced  by  subtracting  the  value  of  from  either  r„(0)  in  the 
autocorrelation  sequence,  or  from  each  element  in  the  principle  diagonal  of  the 
autocorrelation  matrix.  The  difficulty  with  this  method  is  in  determining  the  value  of  the 
noise  power,  r/2.  Additionally,  in  practice  due  to  limited  data,  only  an  estimate  of  the 
autocorrelation  matrix  is  available  in  which  case  equations  (6.11)  to  (6.13)  will  only 
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approximately  be  true.  In  this  case  noise  will  affect  not  only  rii(O),  but  all  other  estimated 
autocorrelation  lags  as  well. 

6.1  EIGEN  ANALYSIS 


Further  insight  into  the  noise  problem  can  be  gained  by  examining  the  structure  of 
the  autocorrelation  matrix  in  more  detail.  For  the  case  of  M  signals,  the  signal 
autocorrelation  matrix  can  be  expressed  as, 


u 

R5  =  (6.14) 

m=  0 

where  Pm  represents  the  signal  power  of  the  incident  signal,  and  the  normalized  signal 
vectors  each  have  the  form  given  by 


Sm  — 


1 


[1, 


(6.15) 


The  parameter  is  the  spatial  frequency  corresponding  to  the  bearing  of  the  signal. 
Since  each  matrix  formed  ffom  the  product  of  the  individual  signal  vectors,  s^Sm^,  has 
rank  1,  the  signal  correlation  matrix,  which  is  the  sum  of  these  matrices,  will  have  rank 
equal  to  the  number  of  signals,  M,  or  have  the  full  rank  p+1  if  M  >  p.  Since  the  signal 
bearings  can  only  be  uniquely  solved  if  Af  <  p  (i.e.  at  least  Af  equations  formed  from  the 
rows  or  columns  of  R  are  are  required  to  uniquely  solve  for  the  A/ unknown  bearings),  then 
only  this  case  is  considered  here. 

One  method  of  examining  the  noise  effects  is  based  on  decomposing  the  signal 
autocorrelation  matrix  into  its  component  eigenvalues  and  eigenvectors.  The  eigenvalues, 
A„  and  corresponding  eigenvectors,  v„  of  a  square  matrix,  M,  have  the  property  that, 

Mvi  =  A, Vi,  (6.16) 


H  H 

where  the  quantity  v,  v,  =1.  If  the  matrix  is  Hermitian  (i.e.  M  =  M),  then  the 
eigenvalues  will  be  re^,  and  if  in  addition,  the  eigenvalues  are  all  distinct,  the  eigenvectors 
will  form  an  orthonormal  set.  The  method  used  to  compute  the  eigenvalues  and 
eigenvectors  are  not  discussed  in  this  report  but  can  be  found  in  reference  [6-1]. 

In  the  case  of  the  signal  autocorrelation  matrix  the  eigenvector  expression 
becomes, 


R5V.  =  AiV.-  (617) 

The  eigenvectors  form  an  orthonormal  basis  set  for  the  signal  vectors,  so  that  the  matrix 
R,  can  be  decomposed  in  terms  of  all  its  eigenvalues  and  eigenvectors  as. 


R,= 


Am^mV  n 


(6.18) 


Tn»  0 

in  which  the  eigenvalues  have  been  ordered  in  decreasing  value  (Ao  >  Ai  >  ...  >  Ap),  and  the 
eigenvectors  associated  with  the  non-zero  eigenvalues  are  orthonormal  (i.e.,  v.t^v^  =  1  if 
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i  =  j,  and  v,f*V;  =  0  if  i  #  j).  The  remaining  eigenvectors  (if  more  than  one),  which  are 
associated  with  the  zero  eigenvalues  can  arbitrarily  chosen  to  satisfy  equation  (6.16)  as 
long  as  the  vector  chosen  is  orthogonal  to  all  the  eigenvectors  associated  with  the  non-zero 
eigenvalues  (but  not  necessarily  those  associated  with  the  zero  eigenvectors).  For 
simplicity,  however,  these  eigenvectors  are  assumed  (throughout  the  rest  of  this  report)  to 
be  chosen  to  satisfy  the  same  orthonormal  condition  as  the  eigenvectors  associated  with  the 
non-zero  eigenvectors. 

Since  the  eigenvectors  associated  with  non-zero  eigenvalues  span  the  same 
subspace  as  the  sign^  vectors,  there  will  only  be  M  non-zero  eigenvalues,  that  is, 

=  Aw*2  =  ...  =  Ap  =  0.  Equation  (6.18)  can  be  rewritten  as, 

u- 1 

Rs=  (6.19) 

m-0 

The  first  M  eigenvectors  are  known  as  the  principal  eigenvectors. 

Similarly  the  decomposition  of  the  autocorrelation  matrix  for  the  case  of  signals  in 
white  noise  can  be  represented  by 


ms  0 


(6.20) 


where  the  eigenvalues  ipt  and  eigenvectors  u,  satisfy  the  relationship 

Ru,  =  (6-21) 

Breaking  the  autocorrelation  matrix  down  into  the  signal  and  noise  autocorrelation 
matrices  (from  equation  (6.13))  then. 


(R,  +  r/2I)Uj  =  ipiU,.  (6.22) 

From  this  expression  it  is  apparent  that  the  autocorrelation  matrix  will  be  full  rank  (i.e. 
p+1  non-zero  eigenvectors)  since  the  noise  autocorrelation  matrix,  r/^I,  has  full  rank. 

The  relationship  between  the  eigenvalues  and  eigenvectors  for  the  autocorrelation 
matrix,  and  those  of  the  signal  autocorrelation  matrix  can  be  highlighted  by  rearranging 
equation  (6.22)  as. 


RjO,  =  (V',  -  ri^W  (6.23) 

This  is,  in  fact,  another  form  of  the  eigenvector  expression  for  the  signal  autocorrelation 
matrix  given  by  equation  (6.17).  Therefore, 

tt,  =  v„  (6.24) 

and 

^,  =  K  +  V^-  (6.25) 
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Using  these  results,  equation  (6.20)  can  be  expressed  in  terms  of  the  signal 
autocorrelation  eigenvalues  and  eigenvectors  as, 


m*  0 


(6.26) 


which  can  be  rewritten  as, 

u- 1  p 

R  =  2^(Am  +  7?2)Vmv“  +  (6.27) 

m  =  0  m  =  U 

This  represents  the  decomposition  of  the  autocorrelation  matrix  into  its  associated 
eigenvalues  and  eigenvectors.  The  first  summation  term  of  equation  (6.27)  contains  the 
eigenvectors  that  span  the  signal  subspace  (i.e.  each  of  these  vectors  is  a  linear  combination 
of  the  signal  vectors),  and  in  this  context  the  first  M  eigenvectors  are  said  to  form  the 
signal  subspace.  The  remaining  p-M+ 1  eigenvectors  are  orthogonal  to  the  signal  subspace 
vectors  and  are  said  to  form  the  noise  subspace. 

In  practice,  since  only  an  est  'mate  of  the  autocorrelation  matrix  is  available,  the 
above  analysis  is  only  approximately  true.  In  this  case,  and  given  no  knowledge  about  the 
noise,  the  left  summation  term  of  equation  (6.27)  represents  the  optimum  reduced  rank 
approximation  in  the  least  squares  sense  of  the  signal  autocorrelation  matrix  [6-2]  (see  also 
Appendix  B).  If  the  noise  is  known  to  be  white  Gaussian  in  nature,  an  even  better  estimate 
of  the  signal  correlation  matrix  can  be  obtained  by  subtracting  an  estimate  of  the  r)^  noise 
contribution  from  the  largest  M  eigenvalues  where  the  estimate  of  rj'^  is  obtained  by  the 
averaging  the  p-M+1  sm^lest  eigenvalues.  Since  removal  of  the  only  affects  the  shape  of 
the  estimated  spectrum,  not  the  location  of  the  peaks,  this  additional  step  is  not  normally 
performed  in  practice. 

6.2  SINGULAR  VALUE  DECOMPOSITION 


Division  of  the  autocorrelation  matrix  into  signal  and  noise  subspaces  can  be 
accomplished  using  eigenanalysis  techniques  as  discussed  in  the  last  section.  A  similar 
division  may  also  be  made  in  the  data  matrix  using  singular  value  decomposition  (SVD) 
techniques.  According  to  th?  SVD  theorem  [6-3],  an  arbitrary  m»n  complex  valued  matrix 
A  of  rank  K  can  be  decomposed  in  terms  of  the  orthonormal  left  singular  vectors, 
uo,  Ui,  tt2,  ...,  Uff-i,  the  orthonormal  right  singular  vectors,  vo,  ▼!,  ▼2,  .  .,  ^k-u  and  the 
positive  real  singular  values,  tro,  (ti,  cr?,  ...,  ax-i,  as, 

K-  1 

A  =  ^<T,Ujv”.  (6.28) 

I  =  0 

where  u,  and  v,  are  m  and  n  element  vectors  respectively,  and  the  singular  values  are 
arranged  in  decreasing  order  (i.e.  cto  >  o’!  >  ...  I  <^*-1  >  0)-  The  relationship  between  the 
singular  vectors  and  the  matrix  A  can  also  be  expressed  as, 


Av,  =  (7,u,  and  A  u,  =  <r,v,. 


(6.29) 


The  relationship  between  SVD  and  eigenanalysis  can  be  •’hown  through  the 
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following  derivation: 


AVi  =  (TjUj 

(Av,)“Avi  =  ((7jUj)Viiii 
H  .  H  .  2  H 

VjA  Avj  =  (Ti  UjUi 

H.H.  2  H  ,  .  H  H 

VjA  Avj  =  cTj  VjVj  (since  UjUj  =  VjVi  =  1) 

A^Avi  =  (Tj  ^Vi .  (6.30) 

Similarly, 

AA^Uj  =  (JiXii.  (6.31) 

From  the  above  two  relationships,  it  is  apparent  that  the  vectors  Vi  are  eigenvectors  of 
the  matrix  A”A  and  the  vectors  u,  are  eigenvectors  of  the  matrix  AA'^.  The  singular 
values  Oi  are  the  positive  square  root  of  the  corresponding  non-zero  eigenvalues  of  either 
matrix. 


If  SVD  is  performed  on  the  signal  data  matrix  then  the  singular  value 
decomposition  of  this  matrix  has  the  form. 


u-  \ 

X5=  2cr,„u,„vi.  (6.32) 

m-O 

For  the  case  of  signals  in  white  noise,  a  similar  decomposition  of  the  data  matrix  can  be 
performed  which  results  in. 


m  =  0 


(6.33) 


where  6^  represents  the  singular  values,  and  w,„  and  represent  the  right  and  left 
singular  vectors  respectively.  However,  since  Vm  is  an  eigenvector  of  the  signal 
autocorrelation  matrix  and  Zm  is  an  eigenvector  of  the  autocorrelation  matrix,  then  from 
the  previous  discussion  on  eigenanalysis  in  Section  6.1  (equation  (6.24)),  Vm  = 

Similarly  it  can  be  shown  that  =  w^.  Therefore 


X  =  (6.34) 

m-O 

Noting  that  am  and  6m  are  the  positive  square  roots  of  the  eigenvalues  of  the 
signal  autocorrelation  matrix  and  autocorrelation  matrix  respectively  (see  equation  (B.20) 
in  Appendix  B),  then  from  equation  (6.25) 

6m=  am+  ■q,  (6.35) 

where  jg  the  noise  power.  Hence  equation  (6.34)  becomes. 
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M-  1 


(6.36) 


X  =  Ufflv”  + 

m=0  m=W 

In  this  form  it  is  apparent  that  the  singular  vectors  associated  with  the  M  largest  singular 
values  of  X  (i.e.  the  vectors  used  in  the  left  summation  term)  are  a  basis  for  the  signal 
subspace,  and  by  default,  th  i  singular  vectors  associated  with  the  p-M+1  smallest 
singular  values  of  X  (i.e.  the  vectors  used  in  the  right  summation  term)  form  the  noise 
subspace. 


In  practice,  since  only  a  limited  amount  of  data  will  be  available,  X“X  will  only 
be  an  estimate  of  the  true  autocorrelation  matrix.  As  a  result,  the  previous  analysis  will 
only  be  approximately  true. 


6.3  QR  FACTORIZATION 

Giv'^n  matrix  A  there  is  a  factorization  such  that 

A  =  QU,  (6.37) 

where  Q  is  an  unitary  matrix  (i.e.  QmQ  =  I)  and  U  is  an  upper  triangular  matrix. 
Eigenanalysis  and  singular  value  decomposition  methods  are  often  performed  based  on 
iterative  QR  approaches  to  solve  for  the  eigenvalues  and  eigenvectors  of  a  matrix. 

An  interesting  alternative  to  the  iterative  approach  is  to  perform  an  approximate 
decomposition  based  on  only  a  single  QR  decomposition  [6-4].  For  example,  performing  QR 
factorization  on  the  signal  autocorrelation  matrix  results  in, 

K,  =  QU,  (6.41) 

which  can  be  expanded  as 

So  =  UooQo  (®-^2) 

8i  =  Moiqo  +  UiiQl 


Sp  —  'aopQo  "h  'JilpQl  -h  ...  4"  ^ppQp) 

where  s,  and  q,  represent  the  columns  of  R,  and  Q  respectively,  and  represents  the 
elements  of  U  .  If  the  signal  matrix  has  rank  M,  then  oidy  the  first  M  columns  of  R, 
will  be  linearly  independent.  That  is,  the  vectors  s^,  Smu  Sm*2)  Sp  will  be  a  finear 
combinatio”  of  the  first  M  columns  of  R,  (so,  Si,  S2,  ...,  Sw-i).  Since  the  vectors 
Qo,  qi,  q2,  qp  are  orthogonal  (by  definition),  then  the  columns  8*f,  Sm+i,  Sm*2,  •■.,  Sp  will  be 
will  ^so  be  a  linear  combination  of  the  vectors  qo,  qi,  q2,  •••,  qiM  only,  and  therefore  the 
values  of  Uy  =  0  for  i  >  M.  From  this  it  is  apparent  that  the  first  M  columns  of  the 
matrix  Q  form  an  orthonormal  basis  set  for  the  signal  autocorrelation  matrix. 

Similarly  the  QR  factorization  of  the  autocorrelation  matrix  results  in, 

R  =  PV,  (6.43) 
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where  R  is  the  autocorrelation  matrix,  P  is  an  unitary  matrix  and  V  is  an  upper 
triangular  matrix.  Again,  this  can  be  expanded  as. 


ro  =  vooPo  (6.44) 

ri  =  uoipo  +  t/iipi 


Tp  =  VOpPo  +  %pl  +  ...  +  VppPp, 

where  ij  represents  the  columns  of  the  autocorrelation  matrix  R.  In  the  case  of  white 
noise  the  relationship  given  in  equation  (6.13)  may  be  used.  Applying  this  to  the  set  of 
equations  (6.44)  then. 


So  +  772(1,  0,  0,  ...,  0]^  = 
si  +  772(0,  1,  0,  ...,  0]’’  = 

7700P0 

VOlPO  +  UllPl 

(6.45) 

8p  +  772(0,  0,  0,  ...,  l]’^  =  ITOpPO  +  VlpPl  +  ...  +  VppPp, 
and  finally  expanding  Sj  in  terms  of  the  relationship  given  in  equation  (6.45)  then 

772(1,  0,  0,  ...,  0]^  +  Uooqo 

T 

772(0,  1,  0,  ...,  0]  +  tioiqo  +  liuqi 

=  VooPO 

=  77oiPo  +  VllPl 

(6.46) 

772(0,  0,  0,  ...,  1]^  +  uopqo  +  uipqi  +  • 

..  +  Up^p  =  VOpPO  +  VipPi  +  ...  +  VppPp. 

By  inspection  the  first  M  basis  vectors  of  the  autocorrelation  matrix  will  not  equal  the 
first  M  basis  vectors  of  the  signal  correlation  matrix,  a  necessary  condition  if  these 
vectors  are  to  be  used  to  create  a  signal  and  noise  subspace.  However,  this  problem  can  be 
overcome  if  the  contribution  of  noise,  is  first  subtracted  from  the  main  iagonal  of  the 
autocorrelation  matrix  before  performing  the  QR  factorization.  In  this  case  the  vectors  q;, 
and  p;  will  be  equal. 

An  estimate  of  the  noise  can  be  made  based  on  the  empirically  obtained  result  that 
I  Vjil  ~  Tf^  for  M  <  t  <  p.  (6.47) 

Rearranging  in  terms  of  and  averaging  over  the  indicated  range  of  the  index  i,  then 

t  -M 

Since  these  results  are  based  on  the  asymptotic  case  where  p-*oo,  the  noise  variance 
estimate  will  be  less  accurate  as  p  decreases. 
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The  ordering  of  the  columns  of  R  is  also  important.  The  optimum  ordering  is  not 
necessarily  the  natural  ordering  resulting  from  the  operation  X^QC,  but  rather  the  ordering 
where  the  first  M  columns  are  selected  to  be  those  which  are  closest  to  being  mutually 
orthogonal.  In  the  case  of  closely  spaced  signals,  the  ordering  of  the  columns  can  be 
determined  in  steps.  At  the  first  step  the  first  column  is  chosen.  In  successive  steps,  the 
column  which  maximizes  the  minimum  distance  (in  terms  of  the  column  indices)  between 
it  and  all  the  other  previously  selected  columns  is  chosen.  The  ordering  of  the  remaining 
columns  is  not  important.  For  example  in  the  case  of  M=2,  the  optimum  ordering  is 
achieved  by  selecting  the  first  column,  then  selecting  the  last  column  since  it  is  farthest 
from  the  first  column.  In  the  case  of  M=3,  the  first  two  colunms  are  selected  the  same  way 
as  for  M=2,  and  the  third  column  will  be  selected  from  a  center  column  (e.g.  for  p=6  the 
ordering  would  be  0,6,2  and  for  p=7  the  ordering  would  be  0,7,3  or  0,7,4). 

The  removal  of  the  noise  power  from  the  main  diagonal  of  the  autocorrelation 
matrix,  followed  by  the  reordering  of  the  columns  of  the  matrix,  are  important  additions  to 
the  QR  method  and  result  in  performance  almost  equal  to  that  of  the  more 
computationally  demanding  eigenanalysis  or  SVD  methods. 

The  noise/signal  subspace  division  is  made  by  considering  that  the  signals  can  be 
completely  represented  by  the  M  column  vectors  of  the  matrix  Q  and  consequently  these 
vectors  form  the  basis  of  the  signal  subspace.  The  remaining  p-M+1  column  vectors  in  Q 
form  the  basis  of  the  noise  subspace. 

In  practice,  when  only  a  limited  amount  of  data  is  available  for  estimation  of  the 
autocorrelation  matrix,  the  preceding  analysis  is  only  approximately  true. 
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7.0  DEALING  WITH  THE  REDUCED  RANK  PROBLEM 


7.1  THE  AUGMENTED  AUTOCORRELATION  MATRIX 

In  Section  5  a  number  of  DF  estimators  were  discussed  (AR,  ME,  LP,  MV,  and 
TN)  which  involve  computing  the  coefficients  of  an  all  pole  filter.  One  common  approach 
to  computing  the  spectrum  of  each  of  these  estimators  involves  the  inversion  of  the 
augmented  normal  autocorrelation  matrix  (R'O-  'I'his  approach  assumes  that  the 
autocorrelation  matrix  R  is  full  rank  and  nonsingular.  In  general  practice  this  will  be  true 
since  uncorrelated  sensor  noise  (usually  modelled  as  complex  white  Gaussian  noise)  will 
ensure  that  the  autocorrelation  matrix  is  full  rank  and  invertible. 

7.1.1  The  General  Solution 

In  Section  6  methods  to  estimate  the  signal  autocorrelation  matrix  R,  .?ere 
discussed.  Using  various  decomposition  methods  to  analyze  the  vector  space  spanned  by 
the  columns  (or  rows)  of  the  autocorrelation  matrix  it  was  shown  that  it  is  possible  to 
divide  this  vector  space  into  a  signal  subspace  and  noise  subspace.  Using  the  matrix  formed 
from  the  signal  subspace  as  the  signal  autocorrelation  matrix  estimate,  the  idea  was  to 
enhance  all  pole  estimator  performance  by  removing  the  noise  before  estimating  the  filter 
coefficients.  The  first  difficulty  with  this  approach,  however,  is  that  the  signal 
autocorrelation  matrix  (and  its  estimate)  has  rank  Af,  so  if  M<p,  this  matrix  is  not 
invertible. 

In  reality  the  difficulty  arises  since  p+1  coefficients  are  being  used  to  determine  M 
signal  bearings.  In  other  words  there  are  more  coefficients  than  necessary  (assuming  M<p) 
and  as  a  result  an  infinite  number  of  solutions  are  possible.  This  difficulty  can  be  overcome 
by  examining  the  all  pole  estimators  in  more  detail.  A  common  feature  of  all  these  methods 
is  the  minimization  of  the  filter  output  variance  given  by  the  equation 

=  aHRa,  (7.1) 

subject  to  some  constraint  (dependent  on  the  estimator  used)  on  the  choice  of  the  filter 
coefficient  vector  a  (the  vector  from  which  the  DF  spectrum  is  generated). 

To  incorporate  the  constraint  into  the  function  to  be  minimized,  the  Lagrange 
multiplier  technique  can  be  used.  Based  on  this  technique  a  new  function  is  defined  such 
that 

F=  aHRa+ 5c(a),  (7.2) 

where  the  function  c(a)  =  0  when  the  constraint  on  a  is  satisfied.  Minimizing  F  with 
respect  to  a^^  results  in  an  equation  of  the  form, 

Ra+&  =  0,  (7.3) 

where  the  constraint  vector  c  is  defined  as 
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Minimizing  F  with  respect  to  6  results  in  the  original  constraint  given  by, 

c(a)  =  0.  (7.5) 

The  solution  of  the  coefficient  vector  a  can  then  be  determined  by  solving  equations  (7.3) 
and  (7.5). 

A  feature  of  the  derivation  of  the  solution  coefficient  vector  a  is  that  it  involves 
the  term  Ra  (see  equation  (7.3)).  If  R  is  invertible,  a  is  easily  isolated  by  premultiplying 
the  term  by  R-^.  In  the  case  where  R  is  not  full  rank,  then  determination  of  a  must  be 
handled  differently. 

The  effect  of  the  constraint  on  the  minimization  of  equation  (7.2)  is  to  reduce  the 
number  of  degrees  of  freedom  in  the  choice  of  the  coefficient  vector  a  by  1.  In  otherwords,  if 
Ris  a  (p-fl)«(p+l)  autocorrelation  matrix,  then  the  choice  of  the  vector  a  will  have  p 
degrees  of  freedom.  The  optimum  choice  for  the  vector  a  is  the  one  wich  results  in  the 
lowest  variance;  ideally  -  q.  However,  if  the  autocorrelation  matrix  is  full  rank,  the  loss 
of  a  degree  of  freedom  means  that  it  is  impossible  to  choose  a  vector  a  which 
simultaneously  satisfies  the  given  constraint  and  sets  the  output  filter  variance  to  zero. 

This  would  require  a  vector  with  p+1  degrees  of  freedom. 

In  the  case  where  the  estimated  signal  autocorrelation  matrix  Rj  is  used  in  place  of 
R,  the  situation  changes  since  R,  is  not  full  rank  (assuming  M<  p)  but  has  rank  M.  Since 
this  reduces  the  number  of  degrees  of  freedom  required  from  p-fl  to  Af,  the  vector  a  can  be 
chosen  to  satisfy  both  the  constraint  and  the  relationship 

a«R^  =  0.  (7.6) 

Expanding  R,  in  terms  of  the  signal  data  matrix  X„  equation  (7.6)  can  be  rewritten  as, 

a%X”a=  |a”X,|^  =  0.  (7.7) 

This  equation  may  only  be  satisfied  if 

'O' 

0 

X»a=  0  .  (7.8) 

6 


64 


Premultiplying  this  result  by  X,  then, 


X,X5a  =  = 


0 

0 

0 


0 


(7.9) 


Assuming  that  M  <  p,  then  there  are  an  infinite  number  of  solutions  to  this 
equation.  The  family  of  all  possible  solutions  may  be  determined  by  examining  the 
eigenstructure  of  R,.  Using  the  notation  developed  in  Section  5.1,  then. 


R,  =  ^A,„v,„y;,  (7.10) 

m=  0 

where  Am  represents  the  eigenvalues  of  R  (and  R,)  ordered  so  that  Ao  >  Ai  >  ...  >  Ap  and  Vm 
represents  the  corresponding  eigenvectors  of  R  (and  R,).  Additionally  R,  has  rank  M  <  p, 
so  that  Xu  =  XuA  =  =  Ap  =  0. 

Since  the  eigenvectors  form  an  orthonormal  basis  set  (vi>^;  =  0  for  t  ^  j,  and 
Vit^Vj  =  1)  and  R,  is  formed  from  the  first  Af  eigenvectors,  then  the  vector  a  must  be 
orthogonal  to  these  M eigenvectors  to  satisfy  equation  (7.9).  Thus  the  vector  a  is  a  linear 
combination  of  the  last  p-M+1  eigenvectors  (v^,  vm*i,  yu*2,  ■■■,  Vp).  One  possible 
representation  is  given  by 


a  =  (  J^mVmVm  )q,  (7.11) 

m-U 

where  the  values  of  6m  are  arbitrarily  chosen,  and  the  vector  q  is  chosen  so  that  a  satisfies 
the  constraint  condition  in  equation  (7.5).  From  equation  (7.11)  a  new  matrix  can  be 
defined  as 


Rl  =  (7.12) 

m-U 

where  the  symbol  T  is  defined  here  as  the  eigeninverse  operator  which  sets  all  the  nonzero 
eigenvalues  of  the  matrix  Rj  to  0  and  all  the  zero  eigenvalues  to  some  arbitrary  value.  The 
choice  of  the  new  eigenvalues  will  be  discussed  later  on  in  Sections  7. 1.2.1  to  7. 1.2.3.  Note 
that  R,i  is  formed  from  the  eigenvectors  of  the  noise  subspace  of  R.  Using  this  new 
representation,  equation  (7.11)  becomes 


a=Rjq.  (7.13) 

The  selection  of  a  suitable  constraint  vector  q  and  subsequent  derivation  of  a  may 
be  simplified  considerably  by  noting  that  when  the  autocorrelation  matrix  has  full  rank, 
equation  (7.3)  can  be  rewritten  as 


a  =  -R'‘6c. 


(7.14) 
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By  setting  q  =  -^c  in  equation  (7.13)  the  solution  of  the  coefficient  vector  a  can  be 
determined  by  solving, 


a  =  -Rl5c,  (7.15) 

subject  to  the  constraint 

c(a)  =  0.  (7.16) 

Solving  these  equations  is  identical  to  the  solution  procedure  for  the  full  rank  case 
(involving  equations  (7.3)  or  (7.14),  and  (7.5)),  except  that  the  matrix  has  been 
exchanged  for  R-i.  In  other  words,  the  noise  removal  techniques  discussed  in  Section  6  can 
be  used  to  enhance  any  of  the  all  pole  DF  estimators  described  previously  in  Section  5 
simply  bv  replacing  R  i  by  the  matrix  RjT.  For  example,  the  enhanced  version  of  the  MV 
method  (see  equation  (5.64))  is  given  by 


5(^)  = 


1 

e“Rje 


(7.17) 


7.1.2  The  Eigeninverse  Solution 

Up  to  this  point  no  criteria  has  been  given  for  the  selection  of  the  eigenvalues  of 
the  eigeninverse  solution  RJ.  The  choice  of  these  eigenvalues  has  an  effect  on  the  size  and 
location  of  spurious  peaks  in  the  DF  spectrum,  and  in  the  case  where  the  estimated 
autocorrelation  matrix  is  used,  it  also  affects  the  location  of  peaks  corresponding  to  actual 
signal  bearings  as  well.  As  a  result,  the  choice  of  the  eigenvalues  can  have  a  significant 
effiect  on  the  accuracy  of  the  DF  estimator. 

Three  different  approaches  are  examined  in  the  next  three  sections. 


7. 1.2.1  The  White  Noise  Approach 

If  the  eigenvalues  of  the  eigeninverse  solution  Rj^  are  chosen,  the  resulting  DF 
spectrum  will  have  M  peaks  corresponding  to  the  bearings  of  the  M  actual  signals,  and  up 
to  jhM  spurious  peaks  with  arbitrary  locations.  A  useful  solution  is  one  that  uniformly 
suppresses  these  spurious  peaks  or  whitens  the  spectrum.  This  reduces  the  chances  of 
confusing  a  spurious  peak  with  an  actual  signal  peak. 


Since  the  eigenvalues  of  a  white  noise  process  are  all  equal,  then  the  estimated 
spectrum  can  be  whitened  by  setting  the  nonzero  eigenvalues  of  R,t  to  1  [7-1].  Defining  this 
in  this  report  as  the  whitened  eigeninverse,  it  is  expressed  in  terms  of  equation  (7.12)  as. 


RI= 


(7.18) 


An  equivalent  definition  using  the  Moore-Penrose  pseudoinverse  [7-2]  is  given  by 


Rl=I-R,R!, 


(7.19) 


where  the  #  is  used  to  represent  the  pseudoinverse  operation.  The  pseudoinverse  of  a 
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matrix  can  be  defined  in  terms  of  its  singular  value  decomposition.  Using  an  arbtrary 
matrix  Q  as  an  example,  then 


m=  0 


(7.20) 


where  K  is  the  rank  of  Q  land  is  less  than  or  equal  to  the  smallest  dimension  of  Q),  am 
represents  the  singular  values  of  Q,  and  and  v^  represent  the  left  and  right  singiilar 
vectors  of  Q  respectively. 

Alternatively,  under  the  special  condition  where  the  signal  data  matrix  has 
dimensions  p*M,  the  whitened  eigeninverse  can  also  be  defined  as. 


rI=  I-X,(X”X,)-‘X!J  (7.21) 

7. 1.2.2  Approach  of  Johnson  and  Degraff 

A  second  approach  proposed  by  Johnson  and  Degraff  [7-3]  to  the  selection  of 
eigenvalues  is  based  on  the  solution  which  results  when  the  autocorrelation  matrix  is 
inverted  first,  and  the  signal/noise  subspace  division  performed  afterwards.  The  result  of 
the  matrix  inversion  is  given  by. 


«- 1 


a‘  = 


1  H 
VmVm 


t 

m-U 


1  H 


(7.22) 


where  the  first  term  represents  the  signal  subspace,  and  the  second  term  represents  the 
noise  subspace.  From  the  analysis  given  previously,  the  enhanced  form  of  R  *  lies  entirely  in 
the  noise  subspace,  in  which  case  the  signal  subspace  eigenvalues  are  set  to  0.  That  is 


R 


T 

s 


H 

m- 


(7.23) 


7.1.2.3  Approach  of  Wax  and  Kailath 

The  third  approach  proposed  by  Wax  and  Kailath  [7-4]  is  based  on  an  examination 
of  the  maximum  likelihood  solution  to  the  bearing  estimation  problem  for  a  single  signal. 
The  maximum  likelihood  method  seeks  to  find  the  signal  which  best  fits  the  data.  For  a 
single  signal,  assuming  zero  mean  Gaussian  noise,  the  maximum  likelihood  solution  reduces 
to  a  least  squares  minimization  which  is  given  by  the  following, 

t  =1 

where  the  minimization  of  is  performed  with  respect  to  c(i)  and  <(>,  the  vector  x(z) 
represents  a  single  sample  of  the  sensor  data  measured  at  the  time  instance  represented  by 
t,  c(t)  represents  the  phase  and  amplitude  of  the  signal  measured  at  sensor  0,  and  e  is  the 
steering  vector  (i.e.  it  determines  the  signal  phase  and  amplitude  at  each  of  the  other 
sensors  as  a  function  of  <p  relative  to  sensor  0). 


(7.24) 
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Minimizing  equation  (7.24)  with  respect  to  c(i)  results  in  the  expression 

=  (7.25) 

e”e 

Replacing  eft)  by  this  expression  and  noting  that  e”e  =  N,  then  the  minimization  of  can 
be  expressed  as  a  function  of  (/>  only.  That  is,  equation  (7.24)  becomes, 

«==  (7-26) 

i  =1 

which  can  be  expanded  as, 

62=  ^x«(t)x(t)  -  (^-;|^)^x“(t)eeMt).  (7.27) 

« =1  1=1 

Given  that  the  above  expression  is  real  valued  and  greater  than  or  equal  to  zero, 
and  since  the  first  term  is  not  a  function  of  0,  then  minimizing  this  expression  is  equivalent 
to  maximizing  the  second  term.  Ignoring  the  constant  multiplier  term,  the  maximum 
likelihood  solution  can  be  found  by  maximizing, 

T 

2]x”(i)ee”x(t).  (7.28) 

« =1 

Additionally  since  e“x(i)  is  a  scalar  value,  the  above  expression  can  be  rewritten  as, 

T 

^e*!x(t)x”(i)e. 

« =  1 

The  covariance  estimate  of  the  autocorrelation  matrix  for  multiple  samples 
defined  by, 

R  =  T 

« 

Using  this  definition  for  R,  the  maximum  likelihood  solution  can  be  determined  by 
maximizing  the  function, 

e«Re, 

which  is  in  fact  the  classical  spectral  estimator.  In  terms  of  the  eigenvalues  and 
eigenvectors  of  the  covariance  estimate  of  R,  this  equation  can  be  reformulated  as, 


(7.29) 
is 

(7.30) 


(7.31) 
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rt  -  1 

e“(  Aovovo  +  ^  AmVmVm  )e. 


(7.32) 


Using  the  fact  that, 


if  -  1  ^ 

^  ^  fmfm  ~  I> 


(7.33) 


and  expressing  vov8  in  terms  of  the  other  eigenvectors,  then  equation  (7.32)  becomes, 

AT-  1 

e”(  Aol  -  2^(Ao  -  A„)vmvII  )e.  (7.34) 


Since  the  term. 


e^  AqIo  =  Aflfi^fe  = 


(7.35) 


is  constant  with  respect  to  its  effect  can  be  ignored.  This  gives, 


-€”(  ^(Ao  -  Am)v,„vl  )e. 


(7.36) 


Finally  maximizing  this  last  expression  is  equivalent  to  minimizing  the  negative  reciprocal 
given  by. 


I 

e"(  ^(Ao  -  A„)v„vll  )e 


(7.37) 


Inspection  of  equation  (7.37)  reveals  that  it  is  in  fact  an  enhanced  form  of  the  MV 
method  for  a  single  signal,  and  that  the  summation  term  in  the  brackets  is  a  form  of  the 
matrix  R,T.  Using  this  form  of  RJ,  Wax  and  Kailath  proposed  an  extension  to  the 
multiple  signal  case  (applying  the  previous  analysis  to  determine  an  exact  solution  for  the 
multiple  signal  case  cannot  be  done  in  any  simple  fashion)  which  is  defined  here  as, 


Rj  —  ^  ”  ABi)V)nVnj, 


(7.38) 


where  Xg  serves  the  same  function  as  Aq  in  equation  (7.38)  and  is  defined  as  the  mean  of  the 
signal  eigenvalues  given  by. 


1 

Ig  =  ^  An,. 


(7.39) 


69 


7.2  THE  NORMAL  AUTOCORRELATION  MATRIX 

Throughout  most  of  this  report,  autoregressive  DF  estimators  based  on  the 
augmented  normal  equation  formulation  of  the  autocorrelation  matrix  have  been  discussed. 
However,  an  alternate  formulation  using  only  the  normal  equations  was  given,  and  is 
repeated  here  as  (see  equation  (5.11)  in  Section  5.2.1), 

RpW  =  -r,  (7-40) 

where  Rp  is  a  pxp  autocorrelation  matrix,  w  is  the  tap  weight  coefficient  vector  (from 
which  the  DF  spectrum  is  derived),  and  r  was  defined  previously  in  equations  (5.13)  and 
(2.5). 


The  signal/ noise  subspace  technique  can  be  used  to  enhance  Rp  in  the  same 
manner  as  in  the  augmented  autocorrelation  case  simply  by  setting  the  smallest  p-M 
eigenvalues  of  Rp  to  zero.  The  result  is  given  by 

u- 1 

R5  =  ^A„V„Vm  ,  (7.41) 

m  =  0 

« 

where  in  this  case  R^  is  the  signal  autocorrelation  matrix  estimate  of  Rp. 

In  equation  (7.40)  the  quantity  of  interest  is  w,  and  it  is  normally  calculated  as 

w  =  -Rp‘r.  (7.42) 

However,  if  Rp  is  replaced  by  R,  this  solution  is  no  longer  applicable  since  R,  is  not  full 
rank  and  therefore  not  invertible.  However,  analysis  of  the  equation 

R,w  =  -r,  (7.43) 

does  yield  some  insight  into  the  solution  to  the  problem.  For  example  w  can  be  represented 
as  a  linear  combination  of  the  eigenvectors  of  the  autocorrelation  matrix  Rp.  One  possible 
representation  is  given  by. 


W=(  ^^,„v,„v„  +  (7.44) 

msO  m-M 

where  in  this  case  q  is  a  column  vector  of  ones  and  the  values  of  6m  are  chosen  so  that  w  is 
a  solution  to  equation  (7.42).  The  above  representation  also  clearly  distinguishes  between 
the  signal  subspace  eigenvectors  (first  summation  term),  and  the  noise  subspace 
eigenvectors  (second  summation  term)  of  the  autc  :orr^ation  matrix. 

Using  this  last  representation  for  w,  and  the  representation  of  R,  given  by  equation 
(7.41),  equation  (7.43)  can  be  expanded  as, 

M-  I  «-  I  p-  1 

RjW  =  (  ^A,„VmVm  )  (  +  ^^rnVmvi  )  q  =  (7  45) 

m-0  m=0  m-M 
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Simplifying,  by  taking  advantage  of  the  fact  that  the  eigenvectors  form  an  orthonormal  set, 
the  result  can  be  expressed  as 


u- 1 

RjW  =  (  2A,„d„,Vmvi )  q  =  -r.  (7.46) 

m  =  0 

From  the  form  of  this  equation,  it  is  clear  that  the  choice  of  values  of  for  M  <  m  <  p  has 
no  effect  on  the  solution  of  equation  f7.43V  In  other  words,  these  values  of  6m  may  be 
arbitrarily  chosen  and  equation  (7.43)  will  still  be  satisfied. 

Based  on  this  observation  it  is  possible  to  derive  the  general  solution  to  equation 
(7.43)  for  w.  The  tap  weight  coefficient  vector  w  can  be  represented  as 

w  =  w,  +  w„,  (7.47) 

where  w,  is  the  part  of  w  which  is  formed  from  the  signal  subspace  of  Rp  in  equation  (7.46) 
and  is  defined  (based  on  equation  (7.44))  as. 


tf- 1 

m3  0 

and  w„  is  the  part  of  w  which  is  formed  from  the  noise  subspace  of  Rp  and ’s  defined  as, 

w„  =  (  )  q-  (7.49) 

m-M 

As  stated  before  the  choice  of  6m  for  w„  is  arbitrary  so  that  equation  (7.49)  represents  the 
general  solution  for  Wn  and  need  not  be  simplified  any  further.  The  choice  of  6m  for  w„  on 
the  other  hand,  is  critical  and  must  be  chosen  so  that  equation  (7.43)  is  satisfied. 

To  solve  for  w„  equation  (7.43)  is  rewritten  as 

R,w  =  R,w,  =  -r.  (7.50) 

Since  R,  is  not  invertible,  a  new  matrix,  derived  from  the  noise  subspace  eigenvectors  (but 
not  eigenvalues)  of  Rp  is  chosen  to  act  as  a  "catalyst”  in  the  solution.  Namdy, 

Rn  =  ^▼mVm-  (7-51) 

This  matrix  is  formed  from  the  noise  subspace  of  R  and  has  the  properties  that, 

Rnw,  =  0,  (7.52) 

and 

R„r  =  0  (7.53) 
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(this  second  property  is  a  result  of  the  fact  that  based  on  equation  (7.46)  r  is  a  linear 
combination  of  the  signal  subspace  eigenvectors  which  are  orthogonal  to  the  eigenvectors  of 
R„).  Using  the  first  property  (equation  (7.52)),  equation  (7.50)  can  be  rewritten  as 

(R,  +  R„)w,  =  -r.  (7.54) 

The  matrix  (R,  +  R„)  is  invertible,  so  the  solution  for  Wj  is  given  by, 


w,  =  -(R,  +  R„)-‘  r.  (7.55) 

Again  using  the  eigenvector  representation,  equation  (7.55)  can  be  expanded  as 


u  ^  ^ 

W,  =  -(  ^  +  ^VmVm  )  T. 


(7.56) 


m=0  m=W 

Recognizing  the  second  summation  term  as  Rn,  and  taking  advantage  of  the  relationship 
expressed  in  equation  (7.53),  then  the  solution  for  w,  simplifies, 

M 


w. 


<  Y  )  '■ 


(7.57) 


Therefore,  the  general  solution  for  w  in  terms  of  the  eigenvalues  and  eigenvectors 
of  the  autocorrelation  matrix  is  given  by, 

u  p-  1 

w  =  -(  ^  ^^mVmVm)q-  (7.58) 

m*  0  m-M 

An  infinite  number  of  solutions  exist  for  w  since  the  values  of  6m  can  be  arbitrarily  chosen. 
7.2.1  The  Pseudoinverse  Solution 


In  the  preceding  discussion  regarding  the  general  solution  of  the  tap  weight 
coefficient  vector  w,  an  infinite  number  of  solutions  were  shown  to  exist.  When  the  true 
autocorrelation  matrix  Rp  is  used,  each  solution  will  result  in  a  DF  spectrum  with  M  peaks 
corresponding  to  the  bearing  of  the  Af  actual  signals,  and  up  to  p-M  spurious  peaks  with 
arbitrary  location.  The  most  useful  solution  is  the  one  that  suppresses  these  peaks. 

In  the  general  solution  of  w  expressed  by  equation  (7.58),  the  first  summation  term 
represents  the  contribution  of  the  signjd  subspace  eigenvectors  or  the  autocorrelation 
matrix  R„  and  remains  invariant  for  all  possible  solutions  of  w.  The  second  term  represents 
the  variable  part  of  w.  From  this  it  can  be  concluded  that  the  si^'-e  and  location  of  the 
signal  peaks  in  the  DF  spectrum  are  controlled  by  the  first  term,  and  the  spurious  peaks 
are  controlled  by  the  second  term.  Consequently  the  size  of  the  spurious  peaks  in  the 
spectrum  can  be  minimized  by  minimizing  the  second  summation  term.  Since  this 
corresponds  to  setting  6m  =  the  desired  form  of  w  is  given  by  [7-5], 

M 

W  =  -(  ]^  )  r  (7.59) 

m*  0 
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The  expression  inside  the  brackets  is  the  pseudoinverse  (previously  defined  by  equation 
(7.20))  of  R,  and  consequently  equation  (7.58)  can  also  be  written  as, 

w  =  -Rf  r,  (7.60) 

where  the  symbol  #  indicates  the  pseudoinverse  operation. 

This  solution  for  w  is  also  called  the  minimum  norm  solution  since  this  solution  of 
w  has  the  minimum  norm  of  all  possible  solutions  where  the  norm  of  w  is  given  by 


|w[  =  ^w“w.  (7-61) 

This  can  be  shown  by  expanding  w  using  the  relationship  expressed  in  equation  (7.58),  and 
squaring  to  get. 


I  W|  ^  ^  j  |^( 

m  =  0  m-M 


M 

0 


1  H  . 

■ymym  )r 


p- 1 

+  ( 

m  =  U 


iVmVn 


)q  (7.62) 


By  taking  advantage  of  the  fact  that  the  eigenvectors  form  an  orthonormal  set,  this 
expression  simplifies  to 


w 


2  =  r'^r 


M 

m*  0 


+ 


q”q 


m  =  U 


(7.63) 


By  inspection  it  is  clear  that  |  wl^  and  correspondingly  the  norm  jwj  are  minimized  when 
6m  =  0.  Therefore  equation  (7.59)  represents  the  minimum  norm  solution  of  w. 
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8.0  ENHANCED  ALL  POLE  ESTIMATORS 

In  the  previous  two  sections  the  techniques  to  remove  a  significant  portion  of  the 
noise  from  the  data  or  autocorrelation  matrix  (Section  6)  and  generate  a  suitable 
replacement  for  the  inverse  autocorrelation  matrix  were  discussed  (Section  7).  Of  particular 
interest  is  the  improvement  that  results  in  all  pole  DF  estimators  when  these  techniques 
are  used.  Enhancing  these  estimators  in  this  manner  is  equivalent  to  using  an  ARMA 
model  without  the  need  to  explicitly  compute  the  MA  filter  coefficients  (since  these 
coefficients  ideally  only  model  the  noise  and  are  therefore  not  required  to  generate  the  DF 
spectrum!.  As  was  discussed  in  Section  4,  the  ARMA  filter  provides  the  best  model  of  the 
sensor  information  in  a  signal  plus  noise  environment. 

MA  estimators  can  also  be  improved  using  the  enhancements  discussed  in  Section 
6.  However,  since  for  small  tactical  sensor  arrays  MA  estimators  are  generally  resolution 
limited  by  the  number  and  spacing  of  sensors,  and  not  the  noise,  these  enhancements  are 
not  discussed  here. 

Enhancements  to  ARMA  estimators  have  been  proposed  [5-1],  but  since  these 
methods  generally  involve  a  computationally  extensive  nonlinear  search  procedure,  they  are 
not  discussed  in  this  report. 

From  the  discussion  in  Section  7  it  was  shown  that  an  infinite  number  of  solutions 
exist  as  possible  replacements  for  the  inverse  autocorrelation  matrix  which  is  a  central 
feature  of  the  all  pole  estimators  discussed  in  Section  5.  In  particular  the  eigeninverse,  the 
weighted  eigeninverse,  and  pseudoinverse  solutions  have  all  been  proposed  as  useful 
solutions.  The  result  is  numerous  possible  enhanced  estimators.  To  limit  discussion  on 
these  enhanced  estimators,  only  those  which  have  appeared  in  the  open  literature  are 
discussed.  These  modified  estimators  are  also  divided  into  two  categories;  linear  prediction, 
and  Capon. 

8.1  ENHANCED  LINEAR  PREDICTION  ESTIMATORS 

The  linear  prediction  estimators  discussed  in  Section  5  are  considered  to  include 
the  Autoregressive,  Maximum  Entropy  and  all  linear  prediction  methods.  In  the  following 
two  sections  two  enhanced  linear  prediction  estimators  are  described.  The  first  is  based  on 
using  the  augmented  normal  equation  formulation,  and  the  second  is  based  on  using  just 
the  normal  equation  formulation. 

8.1.1  Minimum  Norm  Method 

In  the  Minimum  Norm  (MNorm)  method  [8-1]  the  linear  prediction  equation  given 

by, 


Sli{<P)  = 


1 

e^R'uu^R'e 


is  modified  to  become 


,5,VAform(0)  — 


1 

e^Rjuu'^Rje 


(8.1) 


(8.2) 
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where  the  inverse  autocorrelation  matrix  is  replaced  by  the  whitened  eigeninverse  matrix 
(equation  (7.18)).  In  terms  of  the  eigenvalues  and  eigenvectors  of  the  autocorrelation 
matrix,  the  above  expression  can  be  rewritten  as, 

SuSorm{<P)  =  ~y~  •  (^-3) 

e“(  ^  v,„Vm")  uu“  (  ^  v„v„”)e 

m=U  m=U 

The  name  Minimum  Norm  comes  from  the  fact  that  the  filter  coefficient  v'^ctor 
given  by, 

a  =  <r^Iu,  (8.4) 

has  the  minimum  norm  property.  That  is,  although  there  are  an  infinite  number  of 
solutions  for  the  vector  a  which  satisfy  the  equation 

Rja  =  (8.5) 

(where  R,  is  the  estimated  signal  autocorrelation  matrix  defined  in  Section  6),  the  solution 
selected  is  the  one  which  minimizes  the  vector  norm  |  a^aj  subject  to  the  constraint  oo  =  1. 

8.1.2  Modified  Forward-Backward  Linear  Prediction 

Like  the  Minimum  Norm  method  the  Modified  Forward-Backward  Linear 
Prediction  (MFBLP)  method  [8-2]  is  based  on  enhancing  the  forward  backward  linear 
prediction  method  [5-5].  In  the  MFBLP  method,  the  normal  equations  defined  by 

RpW  =  -r,  (8.6) 

are  used,  not  the  augmented  normal  equations.  The  enhanced  solution  for  w  is  given  by, 

w  =  -R#r,  (8.7) 

where  in  this  case  Rf  is  the  pseudoinverse  (see  equation  (7.20))  of  the  signal  autocorrelation 
matrix  derived  from  Rp. 

Once  the  vector  w  has  been  estimated,  the  vector  a  is  determined  using 

a=n.l,  (8.8) 

w 

and  the  DF  spectrum  is  computed  as, 

Smfbli{<P)  =  ^  •  (8.9) 

e”aa^ 

Another  similarity  to  the  Minimum  Norm  method  is  that  the  solution  of  the  vector 
a  determined  by  the  MFBLP  method  has  the  minimum  norm  property.  That  is,  although 
there  are  an  infinite  number  of  solutions  for  the  vector  w  which  satisfy  the  equation 


75 


RjW  =  -r, 


(8.10) 

the  solution  selected  is  the  one  which  minimizes  the  vector  norm  |  w“w| .  If  the  vector  w  is 
a  minimum  norm  solution,  then  so  will  the  vector  a.  From  this  it  would  seem  that  the 
Minimum  Norm  method  and  the  MFBLP  method  are  identical.  However,  the  difference 
between  the  two  methods  lies  in  the  signal/noise  subspace  division.  For  the  MFBLP 
method  this  division  is  performed  on  the  normal  autocorrelation  matrix  Rp  and  for  the 
Minimum  Norm  method  this  division  is  performed  on  the  larger  augmented  autocorrelation 
matrix  R. 

8.2  ENHANCED  CAPON  ESTIMATORS 

The  major  difference  between  the  Linear  Prediction  estimators  and  the  Capon 
estimators  (which  include  the  Minimum  Variance  and  Thermal  Noise  methods)  is  the 
manner  in  which  the  all  pole  filter  coefficients  are  determined  from  the  inverse 
autocorrelation  matrix.  Linear  Prediction  estimators  select  a  single  column  from  this 
matrix  while  Capon  estimators  use  a  linear  combination  of  these  columns. 

Two  enhanced  Capon  estimators  are  described  in  the  following  sections. 

8.2.1  MUSIC  Method 

The  Multiple  Signal  Classification  (MUSIC)  method  [7-lj  can  be  described  as  an 
enhanced  version  of  either  the  Minimum  Variance  or  the  Thermal  Noise  methods.  In  either 
case  the  inverse  autocorrelation  matrix  is  replaced  by  the  whitened  eigeninverse  matrix 
(equation  (7.18)).  The  equivalence  of  the  two  enhancements  (enhanced  Minimum  Variance 
and  enhanced  Thermal  Noise)  is  a  result  of  the  fact  that  the  eigenvalues  of  the  replacement 
matrix  are  set  to  unity  so  that 


Rl  =  (Rl)"'. 

The  resultant  DF  estimator  has  the  form, 


(8.11) 


(8.12) 


In  terms  of  the  eigenvalues  and  eigenvectors  of  the  autocorrelation  matrix,  the  above 
expression  can  be  rewritten  as. 


e”Rle 


Shiusic{<t>)  = 


e  ( 


VmVf] 


)e 


m-M 


(8.13) 


8.2.2  Eigenvector  Method 

The  Eigenvector  (EV)  method  [7-3]  can  also  be  represented  by  the  MUSIC 
estimator  equation  (8.12),  that  is. 


Sev{4>)  — 


1 

e”Rle 


(8.14) 


76 


The  difference  between  the  two  methods  is  in  the  definition  of  RJ.  In  the  EV  method  RJ  is 
the  eigeninverse  matrix  based  on  the  approach  by  Johnson  and  DeGraff  (equation  (7.23)). 
Using  this  definition,  equation  (8.11)  is  no  longer  true,  so  that  the  EV  method  can  be 
viewed  as  an  enhancement  of  the  Minimum  Variance  estimator,  but  not  the  Thermal  Noise 
estimator. 

In  terms  of  the  eigenvalues  and  eigenvectors  of  the  autocorrelation  matrix,  the 
Minimum  Variance  DF  spectrum  can  also  be  defined  as, 

SEv{<t>)  =  - j - ^ - •  (8.15) 

H/  1  H. 

8.3  A  COMPARISON  OF  ENHANCED  DF  ESTIMATORS 

In  this  section,  two  basic  enhanced  estimators  plus  variants  are  compared.  The 
enhanced  estimators  can  be  summarized  as: 

1.  Enhanced  Autoregressive: 

Enhanced  Maximum  Entropy: 

Enhanced  Linear  Prediction: 

2.  Enhanced  Minimum  Variance:  Seikv{<P)  =  — ^ —  (817) 

e“Rle 

Note  that  the  first  category  includes  the  Minimum  Norm  and  Modified  Forward  Backward 
Linear  Prediction  methods  discussed  in  Sections  8.1.1  and  8.1.2  respectively,  and  the 
second  category  includes  the  MUSIC  and  the  Eigenvector  methods  discussed  in  Sections 
8.2.1  and  8.2.2  respectively. 

Variants  include  the  different  approaches  to  determining  the  eigeninverse,  namely, 
the  whitened  eigeninverse,  the  Johnson  and  Degraff  approach,  and  the  Wax  and  Kailath 
approach,  as  well  as  the  approximate  whitened  eigeninverse  solution  based  on  QR 
Factorization.  For  comparison  purposes,  the  whitened  eigeninverse  based  estimators, 
namely  the  Minimum  Norm  method  and  the  MUSIC  method,  are  adopted  here  as  the  basic 
enhanced  estimators. 

The  simulations  that  follow  are  also  based  on  the  identical  data  used  previously  for 
comparing  DF  estimators  in  Section  5.4.  The  optimum  model  order  p  for  spatial  smoothing 
purposes  was  based  on  a  comparison  of  enhanced  estimator  accuracy  versus  the  model 
order,  the  results  of  which  are  shown  in  Figure  8.1.  The  variance  of  the  bearing  errors  for 
each  value  of  p  was  computed  from  300  trials  at  a  fixed  signal  to  noise  ratio  of  50  dB.  With 
the  exception  of  the  variants  using  the  approach  of  Johnson  and  Degraff,  the  model  order 
chosen  was  p  =  5.  For  those  variants  using  the  approach  of  Johnson  and  Degraff,  the 
chosen  model  order  was  p  =  4. 


Selp{<P)  = 


ellJuu^Rle 


(8.16) 
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FIGURE  8.1:  Bearing  error  variance  as  a  function  of  the  model  order  p 


Figure  8.2  compares  the  performance  of  the  Linear  Prediction  method  versus  the 
Minimum  Norm  and  Modified  Forward-Backward  Linear  Prediction  methods,  while  Figure 
8.3  compares  the  performance  of  the  Minimum  Variance  method  with  the  MUSIC  method. 
Threshold  performance  for  all  these  methods  is  dominated  by  the  problem  of  the  peaks  in 
the  estimated  DF  spectrum  corresponding  to  the  40  and  50  degree  bearing  signals  beginning 
to  merge  into  a  single  peak.  In  both  cases  the  enhanced  methods  outperformed  the 
unenhanced  methods  in  terms  of  accuracy  and  threshold  performance.  Interestingly,  in 
comparing  these  results  to  Figure  5.4,  the  Adaptive  Angular  Response  method  still  has  the 
best  threshold  performance  (but  not  accuracy)  of  any  of  the  methods  described  so  far. 

Figure  8.4  compares  directly  the  performance  of  the  Minimum  Norm  method,  and 
the  MUSIC  method.  For  these  simulations,  the  performance  of  the  Minimum  Norm  method 
(and  the  MFBLP  method  based  on  Figure  8.2)  was  slightly  better  than  MUSIC  in  terms  of 
threshold  performance.  Above  threshold  there  was  no  significant  difference  in  the  accuracy 
of  any  of  these  methods.  Other  researchers  have  also  verified  these  results  for  two  signal 
environments  (8-3). 

In  comparing  variants  of  MUSIC  and  the  Minimum  Norm  methods,  involving 
either  the  Johnson  and  Degraff  approach,  or  the  Wax  and  Kailath  approach  to  the 
eigenin verse,  no  significant  differences  were  found  using  the  approach  of  Wax  and  Kailath 
and  the  whitened  eigeninverse  as  shown  in  Figures  8.5(a)  and  (b).  The  approach  of  Johnson 
and  Degraff  was  found  to  have  poorer  accuracy  mainly  diie  to  the  restriction  to  lower 
model  orders  (at  higher  model  orders  accuracy  was  significantly  worse  as  indicated  in 
Figure  8.1).  The  restriction  to  lower  model  orders,  compared  to  the  other  approaches,  also 
limits  the  maximum  number  of  signal  bearings  that  can  be  determined. 

Using  the  QR  factorization  to  compute  an  approximate  whitened  eigeninverse, 
there  was  no  degradation  in  performance  for  the  Minimum  Norm  method  (shown  in  Figure 
8.6  (a)),  and  only  a  slight  degradation  in  threshold  performance  for  the  MUSIC  method 
(shown  in  Figure  8.6  (b)). 
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Signal  to  Noise  Ratio  (dB) 

FIGURE  8.2:  Comparison  of  the  LP,  MNorm,  and  MFBLP  estimators 


Signal  to  Noise  Ratio  (dB) 

PTGURE  8.3:  Comparison  of  the  MV  and  MUSIC  estimators 


FIGURE  8.4:  Comparison  of  the  MUSIC  and  MNorm  estimators 
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FIGURE  8.5:  Comi-arison  of  three  different  eigeninverse  approaches 
(a)  Minimum  Norm,  and  (b)  MUSIC 
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Bearing  Error  Variance  ((leg*dcg)  Beating  Error  Variance  (deg*deg) 


FIGURE  8.6:  Comparison  of  enhanced  estimators  with  their  QR  versions 
(a)  Minimum  Norm,  and  (b)  MUSIC 
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9.0  BEARING  ESTIMATION  USING  THE  ROOT  METHOD 

A  problem  with  the  method  of  searching  the  DF  spectrum  for  signal  peaks  as 
outlined  in  Section  4.4  (called  the  spectral  search  method  here)  is  that  the  steering  vector  e 
constrains  the  solutions  to  the  form  of  a  complex  sinusoid.  In  the  derivation  of  the  DF 
estimators  discussed  in  Section  5,  no  such  constraint  is  placed  on  the  computation  of  the 
coefficient  vector  a,  the  vector  from  which  the  DF  spectrum  is  generated  (see  equation 
(4.29)).  As  a  result,  the  solutions  determined  using  the  spectral  method  are  not  properly 
matched  to  the  vector  a.  Figure  9.1  illustrates  an  example  of  this  where  a  single  peak 
indicates  only  one  signal  although  two  signals  are  present.  A  closer  inspection  of  this  peak, 
does  reveal  a  "kink"  which  indicates  the  presence  of  the  second  signal  but  not  an  accurate 
estimate  of  its  bearing. 


12: - 


Direction  (degrees) 


FIGURE  9.1:  Minimum  Variance  bearing  estimation  using  the  DF  spectrum 
(solid  line)  and  the  root  method  (dashed  lines).  The  actual 
signals  were  at  40  and  60  degrees. 


A  better  approach  is  to  use  a  more  general  steering  vector  defined  as. 


s 


1_ 

i/p+1 


1 

s 

S2 

SP 


(9.1) 


where 


3  = 


(9.2) 


in  place  of  the  steering  vector  Cp  in  equation  (4.29)  (see  also  equation  (2.18)),  that  is. 
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(9.3) 


_ 1_ 

s«aa»s‘ 

The  difficulty  here  is  that  the  search  problem  has  been  changed  from  a  one  dimensional 
problem  in  w,,  to  a  two  dimensional  problem  involving  a;,  and  c.  This  problem  can  be 
simplified  by  noting  that  maximizing  equation  (9.3)  is  equivalent  to  minimizing  the  term 

s”a  =  (9.4) 

n  s  0 

Since  this  is  a  polynomial  equation,  an  equivalent  mathematical  representation  is  given  by, 


s”a  =  do Il(l  - p„s'),  (9.5) 

n  =  i 

where  Pn  represents  the  roots  of  this  equation,  but  will  be  referred  to  as  poles  since  p„  also 
represents  half  the  poles  of  equation  (9.3).  (The  other  half  of  the  poles  of  equation  (9.3)  are 
related  by  the  expression. 


Pn.p  =  -Pn  for  1  <  n  <  p.  (9.6) 

For  the  purposes  of  bearing  estimation,  however,  determination  of  the  poles  given  by 
equation  (9.6)  is  unnecessary.) 

Various  rooting  algorithms  exist,  although  not  discussed  here,  which  are  capable 
of  determining  the  poles  pn  for  1  <  n  <  p  given  the  coefficient  vector  a.  Once  the  poles  have 
been  determined,  the  corresponding  bearings,  based  on  equations  (1.2)  and  (9.2),  and  also 
assuming  the  elevation  angle  ip  =0,  are  given  by 

In  the  previous  example  shown  in  Figure  9.1,  the  location  of  the  two  signals  using 
the  rooting  method  is  shown.  The  improvement  over  the  spectral  method  is  obvious  in  this 
case.  The  improvement  is  also  illustrated  in  the  following  simulation  examples.  The 
generation  of  the  data  is  described  in  Section  5.4  and  is  ^so  the  same  data  as  used  in  the 
examples  in  Section  8.3. 

Figure  9.2  illustrates  the  improvement  of  the  root-MUSIC  and  root-Minimum 
Norm  methods  compared  to  the  spectral  search  version  of  the  Minimum  Norm  method. 

The  main  improvement,  which  is  true  for  any  of  the  all  pole  methods,  is  better  performance 
at  lower  signal  to  noise  ratios  (i.e.  a  lower  threshold).  Interestingly  the  improvement  to  the 
MUSIC  method  using  the  rooting  method  is  greater  than  the  improvement  to  the 
Minimum  Norm  method  so  that  root-Music  outperforms  root-Minimum  Norm.  These 
results  are  consistent  with  observations  made  by  other  researchers  for  the  two  signal  case 
[9-1].  Note  that  the  error  variance  above  the  threshold  is  largely  unaffected. 
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FIGURE  9.2:  A  comparison  of  two  root  methods  and  a  spectral  search  method 


Using  variant  approaches  to  calculate  the  eigeninverse,  the  qualitative  differences 
were  the  same  as  the  spectral  search  case.  That  is,  there  was  no  difference  between  root 
methods  using  the  whitened  eigeninverse  or  the  approach  of  the  Wax  and  Kailath,  whereas 
the  approach  of  Johnson  and  DeGraff  resulted  in  poorer  performance.  Figures  9-3  (a)  and 
(b)  summarize  these  results. 
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FIGURE  9.3:  A  comparison  of  root  based  estimators  using  vairious  eigeninverse 
approaches  (a)  root-Minimum  Norm  (b)  root-MUSIC 
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10  0  DETERMINATION  OF  MODEL  PARAMETERS 


In  all  of  the  DF  algorithms  discussed  so  far,  it  has  been  assumed  that  the  model 
order  parameters  such  as  the  number  of  filter  coefficients,  p  and  q,  and  the  number  of 
signals,  M,  was  already  known.  Generally  this  is  not  true,  so  that  methods  to  estimate 
these  parameters,  which  are  discussed  in  the  following  sections,  are  required.  However, 
before  discussing  these  estimation  methods,  it  is  also  useful  to  explore  the  fundamental 
limits  on  the  maximum  number  of  signals  that  a  DF  array  can  be  used  to  exactly 
determine  the  bearing  of. 

10.1  FUNDAMENTAL  LIMITS  ON  THE  MAXIMUM  NUMBER  OF  SIGNALS 

Previously  the  signal  environment,  in  terms  of  the  sensor  input,  was  defined 
as. 


e 


m=  t 


(10.1) 


In  this  analysis  the  effects  of  noise  are  ignored  since  the  object  is  to  determine  the  absolute 
maximum  number  of  signal  bearings  which  can  be  estimated  as  a  function  of  the  number  of 
sensors. 


At  any  instance  in  time,  each  signal  can  be  completely  described  by  three  real 
valued  parameters,  namely,  amplitude  (represented  by  j  cjt)!),  phase  (a  function  of  both 
Cm{t)  and  Om)  and  direction  (or  correspondingly  the  spatial  frequency  Um)-  Since  these 
parameters  may  be  arbitrarily  chosen  within  a  defined  range,  they  are  completely 
independent.  In  other  words,  if  the  signal  environment  is  to  be  modelled,  the  signal  model 
will  require  a  minimum  of  3Af  real  valued  parameters  to  completely  describe  the 
environment  for  a  single  measurement  sample.  For  multiple  samples  the  number  of 
additional  model  parameters  required  is  dependent  on  the  type  of  signals  (uncorrelated  or 
correlated). 

To  begin  with,  the  case  involving  signals  which  are  all  uncorrelated  in  time  is 
examined.  In  this  case,  the  values  of  0m  and  Um  remain  the  same  from  sample  to  sample 
Isince  they  are  geometry  dependent)  but  the  values  of  the  modulating  envelope  Cm{t) 

(which  is  complex  valued)  will  have  changed.  Assuming  enough  time  has  elapsed  between 
samples  so  that  consecutive  values  of  Cn,(^)  are  independent,  then  each  successive  sample 
will  increase  the  number  of  signal  model  parameters  required  by  2M.  If  T  samples  are  made 
then  the  total  is  given  by, 


3M+ 2M(r- 1).  (10.2) 

In  terms  of  the  input  data  from  which  the  model  parameters  are  estimated,  each 
sensor  input  a:„  can  be  represented  as  a  single  complex  value  or  by  two  real  values,  e.g. 
values  representing  amplitude  and  phase.  For  each  measured  sample,  2 A"  real-valued  input 
parameters  are  available  for  processing.  The  independence  of  these  parameters,  however, 
will  be  a  function  of  the  input  signals,  that  is,  the  number  of  independent  values  cannot 
exceed  the  number  generated  by  the  signal  model.  Mathematically  this  can  be  stated  as, 

N,  =  2NT  forAp>2Ar,  (10.3) 
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otherwise, 


Ni  =  Np  for  Np  <  2NT,  (10.4) 

where  Na  is  the  number  of  independent  real-valued  input  parameters  available  for 
processing;  for  brevity,  these  independent  parameters  will  be  referred  to  simply  as  the 
input  parameters. 

The  condition  =  Np  is  necessary  for  the  solution  of  the  signal  model 
parameters.  This  condition  does  not  necessarily  guarantee  a  unique  solution  since 
ambiguities  may  exist  depending  on  the  signal  model  or  array  geometry  (e.g.  with  linear 
arrays  there  is  a  180  degree  ambiguity).  If  <  Np,  however,  no  unique  solution  exists 
unless  some  form  of  a  priori  knowledge  is  incorporated  into  the  solution  (e.g.  modulation 
type).  Based  on  equations  (10.2)  and  (10.3),  this  condition  is  met  when, 

2NT  >  3M  +  2M{T -  1),  (10.5) 

which  in  terms  of  M  can  be  rewritten  as, 

M<  ^  (10.6) 

N-1 

For  T  >  — this  simplifies  to, 

M<N-1.  (10.7) 

The  case  involving  signals  which  are  all  fully  correlated  in  time  can  be  approached 
by  treating  it  as  a  single  composite  signal.  This  is  motivated  by  the  fact  that  not  only  do 
the  values  of  Om  and  ujm  remain  the  same  from  sample  to  sample  (since  they  are  geometry 
dependent),  but  also  the  ratio  of  the  modulating  envelope  Cm{t)fci{t).  Consequently  only 
one  single  complex  value  (or  two  real  values)  is  required  to  describe  the  change  in  the 
signal  model  parameters  from  sample  to  sample  -  analogous  to  a  single  signal  model  in  the 
uncorrelated  case. 

For  a  single  measurement  sample,  the  composite  signal  will  require  3Af 
real-valued  model  parameters.  Assuming  enough  time  has  elapsed  between  samples  so  that 
consecutive  values  of  Cm{t)  are  uncorrelated,  then  each  successive  sample  will  increase  the 
number  of  real-valued  signal  model  parameters  required  by  2.  If  T  samples  are  made  then 
the  total  is  given  by. 


Np  =  ZM+2{T-l).  (10.8) 

In  terms  of  the  input  data,  the  analysis  is  similar  to  the  uncorrelated  case,  with 
the  exception  that  since  the  composite  signal  changes  only  by  a  complex  multiplier  term 
from  sample  to  sample,  each  sensor  sample  will  be  related  to  the  previous  sample  by  the 
same  constant  multiplier  term.  In  other  words,  after  the  first  sample,  only  2  input 
parameters  are  added  to  the  input  parameter  set  per  sample.  Based  on  this,  the  number  of 
input  parameters  available  after  T  samples  is  given  by, 

Na  =  2N+2{T-l)  forM>^  (10.9) 
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and 


=  3M+ 2(r-l)  ioiM<^. 


(10.10) 


Comparing  the  first  expression,  equation  (10.9),  to  equation  (10.8)  it  is  evident  that  for  the 
stated  conditions  <  Np  regardless  of  the  number  of  samples.  In  comparing  the  second 
expression,  equation  (10.10),  to  equation  (10.8),  the  two  expressions  are  identical,  that  is, 
Nd  =  Np.  Clearly  then,  the  upper  limit  on  the  number  of  correlated  signals  is  given  by, 

(10.11) 

Note  this  is  identical  to  the  uncorrelated  case  where  only  a  single  sample  is  used  ( T  =  1). 

Up  to  this  point  only  the  cases  of  uncorrelated  signals  only,  and  correlated  signals 
only  have  been  considered.  The  more  general  case  consisting  of  a  mixture  of  both  types  of 
signals  can  be  solved  in  a  similar  manner  as  before  by  grouping  the  signals  so  that  each 
group  of  signals  consists  of  signals  (called  fundamental  signals  here)  which  are  correlated 
with  each  other  but  are  uncorrelated  with  signals  from  other  groups.  Each  signal  group  is 
then  treated  as  a  single  composite  signal.  In  the  following  analysis  the  total  number  of 
fundamental  signals  (as  before)  is  represented  by  M,  the  total  number  of  composite  signals 
or  correlated  groups  is  represented  by  K,  and  the  total  number  of  fundament^  signals  in 
each  group  is  M*. 

Following  the  previous  analysis  for  uncorrelated  signals,  the  number  of  signal 
model  parameters  generated  is  given  by, 

Np=ZM+2K{T-l),  (10.12) 

where  the  first  term  and  second  terms  in  this  expression  are  equivalent  to  the  first  and 
second  terms  respectively  in  equation  (10.2). 

The  amount  of  input  data  available  from  the  sensors  will  be  as  given  in  equations 
(10.3)  and  (10.4).  Based  on  the  previous  analysis  of  correlated  signals,  the  condition 

M,<^,  (10.13) 

is  also  imposed  since  adding  more  composite  signals  to  the  problem  will  not  increase  the 
upper  limit  on  Af*. 

The  point  at  which  the  number  of  input  parameters  is  sufficient,  that  is,  =  Np, 
occurs  when 

2NT  >  m  +  2K{T -  1).  (10.14) 

This  expression  can  be  rewritten  in  several  ways  to  reveal  some  of  the  fundamental  limits 
involved.  For  example,  in  terms  of  the  number  of  composite  signals  K,  this  expression  can 
be  written  as 
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(10.15) 


K< 


2NT-m 
5(  T-1)  ’ 


which  for  T  >  1.5M  -  N  +  1  simplifies  to, 


K<  N-1. 


The  maximum  limit  on  M,  based  on  equation  (10.14),  can  be  written  as, 

M< 


(10.16) 


(10.17) 


Inspection  of  this  equation  suggests  no  limit  on  M  if  the  number  of  samples  T  is  increased 
to  infinity.  However,  realizing  that 


K 

(10.18) 

*  =  i 

and  letting  Mk  =  2yV/3  (the  upper  limit  for  Af*)  then  an  upper  limit  on  M  independent  of 
the  number  of  samples  is  given  by. 


M<^^.  (10.19) 

Equations  (10.13),  (10.15),  and  (10.17)  to  (10.19)  are  the  general  expressions  for 
the  mixed  environment  which  place  upper  limits  on  the  values  of  Af*,  K,  and  M 
respectively.  These  limits  reflect  the  fundamental  limitations  placed  on  the  number  of 
signals  that  may  be  estimated  for  any  DF  estimator.  Whether  a  particular  estimator  can 
actually  achieve  these  limits  is  dependent  on  the  modelling  approach  involved.  For 
example,  in  a  signal  environment  where  all  signals  are  uncorrelated  (iiT  =  1)  all  the  DF 
estimators  discussed  in  this  report  are  capable  of  estimating  up  to  N-l  signal  bearings  fthe 
upper  limit)  given  a  sufficient  number  of  samples.  In  the  fully  correlated  environment  \k  = 
Af),  DF  estimators  using  the  covariance  method  to  estimate  the  autocorrelation  matrix  do 
not  achieve  the  theoretical  limit  (only  one  half  this  value)  while  estimators  using  the 
modified  covariance  method  do  (see  Section  3  for  a  discussion  of  these  methods  and  their 
limitations).  None  of  the  estimators  discussed  in  this  report  achieve  the  theoretical  bmits 
in  the  mixed  signal  environment  (1  <  fif  <  Af),  and  are  at  best  limited  to  determining  a 
maximum  of  N-l  signal  bearings. 

10.2  MODEL  ORDER  DETERMINATION 

In  the  following  discussion  only  the  all  pole  methods  discussed  in  Sections  5  and  8 
are  considered  since  MA  methods  have  insufficient  resolution  for  small  tactical  arrays,  and 
ARMA  methods  lead  to  computational  difficulties  which  are  beyond  the  scope  of  this 
report.  This  bmits  the  discussion  to  techniques  which  are  used  to  deternune  the  number  of 
aU  pole  filter  coefficients,  or  filter  order,  p.  Two  types  of  signal  environments  are  also 
considered:  correlated  and  uncorrelated. 

In  the  case  of  uncorrelated  signals  up  to  N-  1  signals  can  be  processed,  assuming  a 
large  number  of  sensor  data  samples  (i.e.  T  »  Af),  without  resorting  to  methods  such  as 
spatial  smoothing.  Since  spatial  smoothing  results  in  a  decrease  in  resolution,  and  tactical 
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arrays  are  generally  small  in  terms  of  numbers  of  antenna,  the  most  desirable  choice  for  the 
filter  order  then  is  p  =  N-  1.  The  exception  to  this  is  for  a  small  number  of  sensor  data 
samples  (i.e.  T  <  ION)  where  the  autocorrelation  matrix  estimates  may  be  unstable  in 
which  case  the  optimum  filter  order  will  lie  somewhere  between  the  maximum  value  of 
N- 1  and  the  optimum  filter  order  for  correlated  signals  discussed  in  the  following 
paragraphs. 

In  the  case  of  correlated  signals,  the  choice  of  filter  order  is  not  nearly  so 
straightforward.  Spatial  smoothing,  which  decreases  the  filter  order,  is  required  if  all 
correlated  sources  are  to  be  resolved.  Decreasing  the  filter  order  also  increases  the  stability 
of  the  autocorrelation  matrix  estimates  since  more  data  is  involved  in  the  computation  of 
each  element  of  the  matrix.  This  can  reduce  the  number  of  spurious  estimates  due  to  noise. 
However,  decreasing  the  filter  order  also  decreases  the  resolution  of  the  estimator. 

Figure  10.1  shows  the  effect  of  changing  the  filter  order  on  the  accuracy  of  the 
bearing  estimates  using  the  Linear  Prediction,  Minimum  Variance,  Thermal  Noise,  and 
•Adaptive  Angular  Response  methods.  The  plots  represent  the  results  of  100  trials  of  a  16 
element  (one  third  wavelength  spacing)  DF  system  used  against  5  coherent  signals  with 
bearings  of  40,  70,  100,  140,  and  160  degrees.  Bearing  errors  at  the  higher  model  orders 
(7  <  p  <  14)  are  a  result  of  the  appearance  of  spurious  peaks  in  the  spectrum. 


~  ~  rooi-Lineai  Prediction 


FIGURE  10.1;  Bearing  error  variance  as  a  function  of  the  model  order  p  for 
basic  DF  estimators 


The  optimum  filter  order  for  the  Linear  Prediction  method  has  been  shown 
analytically  [10-1]  to  be  about  p  =  N/Z  using  dther  the  covariance  or  modified  covariance 
methods  to  estimate  the  autocorrelation  matrix.  Based  on  these  analytical  results  and  the 
simulation  results  shown  in  this  report  (Figures  5.1  and  10.1)  a  choice  for  the  model  order 
of  N/Z  <  p  <  N/2  is  appropriate.  At  higW  filter  orders,  the  extra  coefficients  give  rise  to 
spurious  estimates  in  the  spectrum  which  degrade  the  performance  of  the  estimator.  This 
choice  also  limits  the  number  of  signals  that  can  be  handled  to  Af  <  N/2. 
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Given  the  similarity  between  the  effect  of  model  order  on  the  accuracy  of  the 
Linear  Prediction  method  and  on  the  Thermal  Noise  method  (see  Figure  10.1),  the  choice 
for  the  model  order  is  the  same,  namely,  N/3  <  p  <  iV/2.  The  Minimum  Variance  method 
has  greater  immunity  to  spurious  estimates  so  that  slightly  higher  model  orders  are 
possible  except  that  the  accuracy  of  the  method  decreases  at  these  higher  orders  so  that  a 
sdection  for  the  model  order  of  N/3  <  p  <  Nl2  is  again  appropriate. 

The  Adaptive  Angular  Response  method  is  the  most  sensitive  to  model  order  with 
the  best  choice  being  p=  Moi  slightly  greater.  Where  the  number  of  signals  is  unknown, 
algorithms,  such  as  the  ones  discussed  in  Section  10.3  are  required.  This  significantly 
increases  the  computational  work  load  so  that  the  advantage  of  using  this  method 
compared  to  the  enhanced  methods  may  be  lost. 

In  the  case  of  enhanced  DF  methods,  taking  advantage  of  the  noise/signal 
subspace  division  effectively  eliminates  the  spurious  estimate  problem.  Figure  10.2 
illustrates  the  improvement  of  the  enhanced  estimators  using  the  same  data  as  Figure  10.1. 
Researchers  have  determined  that  for  these  methods  (using  the  whitened  eigeninverse  or 
the  approach  of  Wax  and  Kailath)  the  optimum  filter  order  lies  somewhere  between  0.6iV 
and  O.SiV  [7-5], [10-2].  Based  on  maximizing  the  number  of  signal  bearings  which  can  be 
estimated,  the  best  choice  is  p  =  2A/3.  For  the  methods  using  the  approach  of  Johnson  and 
DeGraff  (e.g.  the  Eigenvector  method),  a  lower  model  order  is  required  due  to  a  decrease  in 
accuracy  at  higher  model  orders  (see  the  example  shown  in  Figure  10.3  at  a  model  order  of 
p  =  10).  This  decrease  in  accuracy  restricts  the  approach  of  Johnson  and  DeGraff  to  lower 
model  orders,  i.e.  p  =  N/2. 


lO* 


FIGURE  10.2:  Bearing  error  variance  as  a  function  of  the  model  order  p  for 
enhanced  DF  estimators 


It  should  be  noted  that  the  optimum  model  orders  described  above  have  been 
considered  for  the  case  where  only  a  single  sensor  sample  is  used  for  bearing  estimation, 
and/or  the  signals  are  fully  correlated.  For  a  large  number  of  samples  {T»  N)  and 
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uncorrelated  signals,  the  optimum  choice  for  the  model  order  is  p  =  ^-  1.  In  land  tactical 
VHF,  correlated  signals  (multipath)  are  a  fact  of  life,  so  that  the  previous  analysis  on 
model  order  selection  is  applicable. 


FIGURE  10.3:  Bearing  error  variance  as  a  function  of  the  model  order  p  for  the 
root-MUSIC  estimators  using  different  eigeninverse  approaches 


10  3  SIGNAL  NUMBER  ESTIMATION 

Enhanced  all  pole  DF  methods  such  as  those  discussed  in  Section  8  require 
knowledge  of  not  only  the  filter  order  p,  but  also  the  number  of  signals  M.  Choosing  M  to 
be  too  small  results  in  poor  resolution  and  choosing  M  too  large  results  in  spurious 
estimates. 

Kailath  and  Wax  have  reformulated  the  Akaike  Information  Criterion  (AIC)  and 
the  Minimum  Description  Length  (MDL)  algorithms  which  were  originally  formulated  for 
estimators  using  the  autocorrelation  method  of  autocorrelation  matrix  estimation 
[10-3], [10-4].  The  reformulated  versions  use  the  eigenvalues  of  the  normal  (not  augmented) 
autocorrelation  matrix.  The  details  of  the  derivation  of  these  algorithms  are  found  in 
reference  [10-5],  and  the  final  results  are  given  by, 


AlC(m)  =  N  (p  -  m)  In 


-I-  m  (2p  -  m), 


(10.20) 


and. 
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MDL(m)  =  N  (p-  m)  In 


t  I 


«  =  mTl 


+  Y  fn  {2p  -  m)  \n(N).  (10.21) 


The  best  estimate  of  the  number  of  signals  M  is  the  value  of  the  parameter  m  which 
minimizes  either  one  of  the  above  functions. 


Studies  have  shown  [10-6]  that  for  a  limited  number  of  data  samples  the  AIC 
algorithm  performs  best.  However,  this  algorithm  is  not  consistent,  that  is,  as  the  number 
of  samples  is  increased  the  probability  of  error  does  not  decrease  to  zero.  The  MDL 
algorithm  is  consistent  so  that  for  large  numbers  of  data  samples  it  performs  better  than 
the  AIC  algorithm. 
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11.0  SUMMARY  AND  CONCLUSIONS 


In  this  report,  a  number  of  the  more  popular  superresolution  DF  methods  were 
examined.  In  particular  those  methods  which  were  based  on  adaptive  filter  models.  Three 
filter  models  have  evolved,  namely,  the  moving  average  filter,  the  autoregressive  filter,  and 
the  autoregressive  moving  average  filter.  In  terupgjif  accuracy  and  computational  speed, 
methods  based  on  the  all  pole  or  autoregressive-nl^  model  have  shown  the  greatest 
promise  and  for  this  reason  were  the  main  focus  of  this  report. 

It  was  shown  that  the  autoregressive  filter  based  methods  can  be  generalized  as  a 
five  step  procedure.  These  steps  include: 

1.  Estimation  of  the  autocorrelation  matrix. 

2.  Division  of  the  autocorrelation  .matrix  into  a  signal  ana  noise  subspace. 

3.  Generation  of  an  enhanced  inverse  autocorrelation  matrix. 

4.  Estimation  of  the  all  pole  filter  coefficients. 

5.  Estimation  of  the  signal  bearings  from  the  filter  coefficients. 

Each  of  these  steps  has  been  the  focus  of  research  on  ways  of  optimizing  the  DF  estimation 
process. 


In  terms  of  the  resultant  DF  accuracy  and  computational  simplicity,  the  modified 
covariance  method  has  been  found  to  be  the  best  choice.  Other  methods  including  the 
autocorrelation  and  the  covariance  methods  result  in  poorer  accuracy,  and  maximum 
likelihood  based  methods,  which  potentially  give  higher  accuracy,  are  currently  too 
computationally  intensive  for  realtime  applications. 

Enhancements  to  DF  estimators  based  on  the  division  of  the  autocorrelation 
matrix  into  a  noise  and  signal  subspace  has  proved  to  be  extremely  useful  for  a  number  of 
reasons.  One  is  the  suppression  of  spurious  l^arings  estimates  which  can  seriously  degrade 
accuracy.  A  second  reason  is  the  ability  to  operate  using  higher  order  filters  which  leads  to 
a  slight  improvement  in  accuracy  and  the  ability  to  determine  a  greater  number  of  signal 
bearings.  A  third  reason  is  that  methods  which  have  been  enhanced  perform  better  at  lower 
signal  to  noise  ratios. 

Two  methods,  eigen  decomposition  and  singular  value  decomposition,  are 
typically  used  to  perform  the  signal/noise  subspace  division.  A  new  autocorrelation  matrix 
estimate  is  then  generated  using  only  the  signal  subspace.  From  a  theoretical  point  of  view, 
the  results  are  identical  using  either  decomposition  method.  An  approximate  method  called 
the  QR  factorization  method  has  also  been  proposed  which  achieves  nearly  the  same 
performance  as  ither  eigen  decomposition  and  singular  value  decomposition,  but  with 
significantly  less  processing  requirements.  This  suggests  that  more  research  is  required  to 
determine  how  accurately  the  signal/noise  subspace  division  must  be  done,  and  whether 
there  are  even  faster  methods  that  could  be  used. 

The  difficulty  with  using  the  signal  subspace  estimate  for  the  autocorrelation 
estimate  is  that  the  resultant  matrix  is  noninvertible.  Since  the  inverse  autocorrelation 
matrix  is  a  central  part  of  all  the  autoregressive  filter  based  estimators,  generation  of  an 
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enhanced  inverse  matrix  is  required.  Four  approaches  have  been  proposed  to  do  this, 
namely  the  pseudo-inverse,  whitened  eigeninverse,  the  approach  by  Johnson  and  Degraff, 
and  the  approach  by  Kailath  and  Wax.  With  the  exception  of  the  approach  by  Johnson  and 
Degraff  the  differences  between  the  various  approaches  were  insignificant.  The  approach  by 
Johnson  and  Degraff  was  found  to  be  restricted  to  smaller  model  orders  than  the  other 
approaches  with  a  resultant  degradation  in  accuracy  and  reduction  in  the  number  of  signal 
bearings  that  can  be  estimated.  From  the  standpoint  of  simplicity,  the  whitened 
eigeninverse  is  the  preferred  approach. 

The  methods  for  determining  the  autoregressive  filter  coefficients  can  be  divided 
into  two  groups:  Capon  estimators  (Minimum  Variance  and  Thermal  Noise)  and  linear 
prediction  estimators  (Autoregressive,  Linear  Prediction,  and  Maximum  Entropy). 

Without  the  signal/noise  subspace  enhancements  (steps  2  and  3)  the  linear  prediction 
estimators  are  characterized  by  higher  resolution  abilities  while  the  Capon  estimators  (at 
least  in  the  case  of  the  Minimum  Variance  estimator)  are  characterized  by  better 
suppression  of  spurious  estimates.  A  third  type  of  estimator,  based  on  a  pole-zero  filter 
model,  had  the  best  threshold  performance  in  the  simulations  performed  for  this  report. 

Of  the  enhanced  estimators,  the  enhanced  linear  prediction  estimators  (e.g. 
Minimum  Norm  and  MFBLP)  have  better  threshold  performance  and  correspondingly 
better  resolution  (when  the  bearings  are  determined  from  the  DF  spectrum)  than  the 
enhanced  modified  Capon  estimators  (i.e.  MUSIC). 

Two  methods  are  currently  in  use  for  estimating  the  bearings  based  on  the  filter 
coefficients.  The  first  and  more  commonly  used  method  is  based  on  a  search  of  the  DF 
spectrum  which  is  generated  based  on  the  transfer  function  of  the  filter.  This  method 
suffers  from  a  loss  in  resolution  for  two  closely  spaced  signals  (in  bearing)  at  lower  signal  to 
noise  ratios.  A  more  recently  popularized  method  is  based  on  determining  the  bearings 
directly  from  the  roots  of  the  polynomial  equation  formed  using  the  filter  coefficients.  The 
result  is  no  merging  of  signal  bearing  as  in  the  spectral  search  case,  and  a  substantial 
improvement  in  performance  at  low  signal  to  noise  ratios.  At  high  signal  to  noise  ratios, 
the  performance  of  both  methods  is  identical.  Based  on  this,  the  root  method  is  clearly 
superior.  Interestingly  the  root  method  brings  the  greatest  improvement  to  the  Capon 
estimators,  that  overall,  the  root-MUSIC  method  had  the  best  performance. 

Other  considerations  affecting  DF  include  estimation  of  the  number  of  signals 
and  the  optimum  number  of  filter  coefficients.  Estimation  of  the  number  of  signals  is 
important  for  the  enhanced  methods  since  underestimating  the  number  of  signals  leads  to  a 
serious  degradation  in  accuracy  and  overestimating  increases  the  possibility  of  spurious 
estimates  which  can  also  seriously  degrade  accuracy.  Two  algorithms,  namely,  Akaike 
Information  Criteria  (AIC)  and  Minimum  Description  Length  (MDL),  have  been  proposed. 
The  AIC  algorithm  has  better  performance  for  a  limited  number  of  sensors  and  sensor  data 
samples.  For  a  larger  number  of  samples  the  MDL  algorithm  performs  better. 

The  choice  of  model  order  is  important  since  this  affects  the  accuracy  of  the 
resultant  bearing  estimates  and  the  number  of  bearings  that  can  be  estimated.  Additionally 
the  basic  DF  methods  are  more  susceptible  to  spurious  estimates  at  higher  model  orders. 

For  the  basic  methods,  the  optimum  model  order  is  generally  the  highest  model  order  for 
which  spurious  estimates  do  not  occur.  Since  this  is  a  function  of  the  number  of  signals 
which  will  generally  be  unknown  (the  AIC  and  MDL  algorithms  require  eigendecomposition 
of  the  autocorrelation  matrix  and  are  therefore  more  appropriate  for  the  enhanced 
methods)  researchers  have  determined  the  best  choice  for  the  model  order  is  between 
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N/Z  <  p  <  N/2  when  the  number  of  signals  is  small  (i.e.  M<  p).  For  methods  which  are 
sensitive  to  model  order,  such  as  the  Adaptive  Angular  Response,  this  criteria  may  not  be 
sufficient  to  guarantee  good  results.  For  the  enhanced  methods  researchers  have  found  a 
choice  of  0.6N  <  p  <  O.SNto  be  appropriate. 

Based  on  the  results  of  comparisons  between  various  methods  in  this  report,  the 
root-MUSIC  estimator  using  the  modified  covariance  autocorrelation  matrix  estimate 
represents  the  best  approach.  Approximate  decomposition  techniques  such  as  QR 
factorization  show  great  promise  in  decreasing  the  computational  requirements  of  the  DF 
estimation  process  without  sacrificing  accuracy,  and  this  avenue  of  research  should  be 
explored  more  fully. 

The  concept  of  signal/noise  subspace  division  has  proved  to  be  a  major  step 
forward  in  DF  estimation.  However  the  original  concept  is  based  on  approximations  which 
are  valid  for  a  large  number  of  sensors.  More  research  in  the  area  of  limited  numbers  of 
sensors  is  necessary. 

It  should  be  recognized  that  the  comparisons  in  this  report  have  focused  primarily 
on  artificial  multipath  type  environments.  The  purpose  of  this  is  due  to  the  fact  that  for 
tactical  DF  systems,  the  problems  introduced  by  multipath  in  the  operational  environment 
are  an  order  of  magnitude  more  serious  than  any  other  error  mechanism  and  remain  largely 
unsolved.  The  word  "artificial"  is  appropriate  since  there  is  very  little  documented  results 
available  on  multipath  measurements  at  VHF  and  UHF  appropriate  for  DF  research  so 
that  in  this  sense  the  multipath  environments  were  contrived.  It  is  quite  likely  that  the 
numbers  of  multipath  signals  in  the  real  environment,  even  for  a  single  transmitter,  will 
easily  exceed  the  number  of  sensors  of  a  tactical  array.  Under  these  conditions,  the 
assumptions  made  in  deriving  the  various  DF  methods  will  be  in  error  and  will  require 
modifications.  However,  many  of  the  techniques  described  in  this  report  will  still  be 
applicable. 

Two  other  somewhat  artificial  assumptions  are  the  assumptions  of  white  Gaussian 
noise  statistics  and  perfect  calibration  of  the  array  sensors.  In  practical  systems  these 
assumptions  are  rarely  true  and  the  result  can  be  a  severe  degradation  in  the  performance 
of  a  DF  estimator.  Although  coloured  noise  can  be  handled  by  prewhitening  the  input 
sensor  data  using  a  preprocessor,  this  requires  knowledge  of  the  noise  statistics  beforehand 
(i.e.  the  noise  only  autocorrelation  matrix).  Determining  unknown  noise  statistics  and 
sensor  calibration  are  both  active  areas  of  research. 

It  should  also  be  recognized  that  the  discussion  in  this  report  has  been  limited  to 
linear  arrays  which  have  specid  mathematical  advantages  compared  to  nonuniform  or 
nonlinear  arrays.  As  a  result  some  of  the  techniques  described  in  this  report  are  only 
applicable  to  uniform  linear  arrays.  This  includes  spatial  smoothing  and  rooting  techniques. 
Additionally  the  modified  covariance  method  is  also  restricted  to  symmetric  arrays.  A 
solution  to  this  problem  has  been  suggested  in  a  paper  by  Friedlander  [11-1]  which 
describes  how  an  arbitrary  array  can  be  interpolated  to  a  linear  array  before  processing. 

The  result  is  that  the  techniques  discussed  in  this  report  can  be  extended  to  arbitrary 
arrays,  albeit  with  some  extra  computational  complexity. 
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13.0  GLOSSARY 


13.1  MATHEMATICAL  CONVENTIONS 

M  a  boldface  upper  case  letter  represents  a  matrix 

V  a  boldface  lower  case  letter  represents  a  column  vector 

s  or  S  a  lower  case  or  upper  case  italic  letter  represents  a  scalar  quantity 

*  complex  conjugate 

T  matrix  (or  vector)  transpose 

H  complex  matrix  (or  vector)  transpose 

estimate 

T  matrix  eigeninverse  (defined  in  section  7.1.2) 

*  matrix  pseudoinverse  (defined  in  section  7. 1.2.1) 

E{}  expectation  of 

real{}  real  part  of 

imag{}  imaginary  part  of 


13.2  SYMBOL  DEFINITIONS 

d  spacing  between  sensors 

A  signal  wavelength 

noise  variance 

T  number  of  time  samples 

M  number  of  signals 

iV  number  of  sensors 

u  spatirJ  frequency  (as  opposed  to  the  usual  definition  for  temporal 

frequency) 

<p  bearing  angle 

(pm  bearing  angle  of  signal  m 

Xm{t)  complex  baseband  output  of  sensor  m  at  time  instance  t 
Xm  same  as  Xm{t)  except  the  time  instance  is  unspecified 

X  sensor  data  vector  (defined  in  section  3) 

X  sensor  data  matrix  (defined  in  section  3) 

rxJ^m)  autocorrelation  lag  m  (defined  in  section  2.1.1) 

R  augmented  autocorrelation  matrix  (defined  in  section  5.2.1) 

Rp  normal  autocorrelation  matrix  (defined  in  section  5.2.1) 

R5  autocorrelation  matrix  formed  mom  signals  only 

Rn  autocorrelation  matrix  formed  from  noise  only 

r,;  element  of  the  autocorrelation  matrix  occupying  the  row  and 

column 

(Xi  ith  singular  value  of  the  data  matrix  ordered  so  that  a*  >  am 

A,  jth  eigenvalue  of  the  autocorrelation  matrix  ordered  so  that  Ai  ^  Ai*i 

Vj  eigenvector  corresponding  to  Aj  or  right  singular  vector  corresponding 

to  ai 

Uj  left  singular  vector  corresponding  to  a* 

e  steering  vector  (defined  by  equation  2.18) 

a  autoregressive  filter  coefficient  vector 

b  moving  average  filter  coefficient  vector 

p  autoregressive  (all  pole)  filter  order  -  p+1  is  also  the  spatial  smoothing 

subarray  size 
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q  moving  average  (all  zero)  filter  order 

5(0)  direction  finding  spectrum 


13.3 


ACRONYMS 

AAR 

adaptive  angular  response 

AIC 

Akaike  information  criterion 

AR 

autoregressive 

ARMA 

autoregressive  moving  average 

BART 

Bartlett 

DF 

direction  finding 

EV 

eigenvector 

MDL 

minimum  description  length 

MUSIC 

multiple  signal  dassification 

MV 

minimum  variance 

TN 

thermal  noise 
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APPENDIX  A  -  POLYNOMIAL  FACTORIZATION 


In  classical  and  superresolution  DF  estimators  functions  of  the  form 


G{uj)  =  e“De, 

are  often  encountered  where  e  is  a  q+1  element  steering  vector  of  the  form, 


e  —  [  1,  c  ,  e  ,  ,  c  j 


-juqd.T 


(Al) 


(A2) 


and  D  is  a  (g+1)  »  (9+I)  Hermitian  symmetric  matrix  (such  as  the  autocorrelation  matrix 
or  its  inverse).  In  expanded  form,  the  function  G{u)  may  also  be  represented  as. 


G{uj)  =  ^  (A3) 

n=-  q 

where 

q  *  n 

diitn,  (■A.4) 

I  *0 

and  is  the  element  located  in  the  row  and  j+n'^^  column  of  the  matrix  D.  Given  the 
properties  of  the  matrix  D,  the  coefficients  On  will  also  have  the  property  that, 

dn  —  tt-n-  (■^^) 

Since  equation  (A3)  represents  a  polynomial  expression  with  an  odd  number  of 
coefficients,  it  can  be  decomposed  into  the  product  of  simpler  polynomial  terms  given  by, 


G{u})  =  rx,(0)  +  1  -  c„c'^). 

n  *  1 


(A6) 


The  polynomial  term  on  the  right  hand  side  of  this  expression  can  be  expressed  in  terms  of 
its  roots  as. 


•  *i(Jd 


+ 1- 


-}ud.  -e 


*  ;ud 


Pin  +  P2n  ■ 


(A7) 


From  inspection  of  equation  (A7),  if  p„  is  a  ’•oot,  then  1/p^  must  also  be  a  root. 
This  fact  can  be  used  to  advantage  although  two  separate  cases  must  be  considered.  In  the 
first  case,  where  |pn|  1)  the  roots  exist  in  pairs  which  must  satisfy  the  inverse  conjugate 
relationsUp.  Under  these  conditions  the  coefficients  c„  are  chosen  so  that 


Pin  =  Pn 


(A8) 
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and 


P2n  =  1/P«*  (A9) 

Using  these  new  relationships  equation  (A7)  becomes, 

-  jud  .  *  j  u  d 

+  1  -  =  (l  -  -i_)  - 

P«  PnPn  +  1 

=  - .  (AlO) 

PnPn  +  1 

In  the  second  case  where  |  p„|  =  1,  the  root  pairs  do  not  necessarily  exist  in  pairs 
since  each  root  of  unit  magnitude  is  its  own  conjugate  inverse,  that  is, 

p,  =  4-  (All) 

Pn 

However,  assuming  that  the  function  G{u)  has  the  property  that 

G(w)  >  0  for  all  w,  •  (A12) 

(such  as  when  G{u)  represents  the  true  spectrum  or  its  inverse)  and  ignoring  cases  where 
different  quadratic  polynomial  terms  have  the  same  roots  (in  that  case  the  c„  are  chosen  so 
that  Pin  =  p2n  so  the  analysis  used  when  |p„|  ^  I  still  applies)  then  it  also  follows  that 


«  *]Ud 

-c„e 


+ 


1  -  c„g  >  0, 


(A13) 


By  inspection  (e.g.  letting  w,  =  0), 


Cn  <  (A14) 

Equation  (A7)  is  a  quadratic  polynomial  so  that  the  two  poles  pi„  and  p2n  can  be 
determinea  from  Cn  using  the  quadratic  equation.  This  gives, 


Pn  = 


(A15) 


Since  |p„|  =  1,  then 


1  _  1  ±  V  1  ~  4CnCn 

2c,  ■ 

Rearranging  this  result  in  terms  of  c„  and  using  the  fact  that  |  c„|  <  1/2,  then 


(A16) 


yj  4CnCn  —  1  *  1  ~  4CnCn 


(A17) 
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Squaring  both  sides, 


4CnCn  —  1  *  2^  1  -  4CnCn  4"  (1  “  4CnCn), 

(A18) 

and  simplifying. 

-(1  -  4c„c^)  =  1  -  4c„4 

(AlS) 

Again  squaring  both  sides 

(1  -  4CnCn)  ~  1  -  4CnCn, 

(A20) 

and  simplifying 

4c„c;i(l  -  4c„c^)  =  0. 

(A21) 

The  solution  c„Cn  =  0  is  not  a  valid  solution  to  equation  (A15)  so 

.  1 

CnC„  =  j. 

(A22) 

Based  on  this  last  result,  equation  (AlS)  simplifies  to. 

(A23) 

Since  there  is  only  one  solution  for  both  pt„  and  p2n,  they  both  must  be  equal.  That  is, 

Pin  —  P2n  —  Pn 

(A24) 

These  roots  also  satisfy  the  relationships  expressed  in  equations  (A8)  and  (A9),  and  as  a 
result  equation  (A  10)  is  still  valid  for  this  case. 

From  the  preceding  analysis,  if  the  condition  >  0  is  met  for  all  u,  and  using 
the  result  given  in  equation  (AlO),  then  equation  (A6)  can  be  factored  into  two  conjugate 
parts  as, 


G{u)  =  r„(0) 


rt 


(1 


"  '  ‘  /  PnPn  +  1 


" “ ‘  V  PnPn  +  1 


Finally,  converting  each  multiplier  term  back  to  polynomial  form,  the  result  is 


(A25) 


0(0,)  =  ^  4.  J  4; 


(A26) 


where  the  coefficients  represented  by  are  a  result  of  combining  terms  with  the  same 
power  of  e. 
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APPENDIX  B  -  OPTIMUM  REDUCED  RANK  MATRICES 


The  problem  to  be  examined  in  this  discussion  is  that  given  an  m  *  n  matrix  A 
of  rank  K,  what  is  the  best  reduced  rank  fless  than  K)  matrix  equivalent  of  A?  The 
definition  of  best  reduced  rank  matrix  can  oe  given  a  mathematiczd  formulation  by 
considering  the  error  matrix, 

E  =  A-B,  (Bl) 

where  B  is  an  m  «  n  reduced  rank  equivalent  of  matrix  A.  The  best  matrix  B,  in  the 
least  squares  sense  minimizes  the  error  power  given  by, 


m-  1  n-  1 

^  ^  CijCi;,  (B2) 

«  =  0  ;  =  0 

where  e,j  represents  an  element  of  the  matrix  E  in  equation  (Bl). 

One  approach  to  solving  this  problem  is  to  consider  an  orthonormal  vector  basis 
set,  Vo,  vi,  ...,  Vif-i  (where  v,  is  an  n  element  column  vector),  used  to  represent  the  rows 
of  the  matrix  A.  For  example,  if  is  an  n  element  row  vector  representing  the 
row  of  the  matrix,  then 


(B3) 

Multiplying  both  sides  of  this  expression  by  V;  then 

ao“v;  =  70;,  (B4) 

H 

ai  V;  =  71;, 
a2  V;  =  72;, 

H 

an»-l  V;  —  7to-1;, 

since  by  definition  v,^*V;  =1  if  i  =  j  and  Vi“V;  =  0  if  j. 

Based  on  the  set  of  equations  (B4)  a  second  vector  set  can  be  defined  such  that, 

Av;  =  (TjU;  (B5) 

where  ui  is  an  m  element  normalized  column  vector,  and  aj  is  a  positive  real  scalar. 

The  matrix  A  can  now  be  decomposed  in  terms  of  the  vectors  Uj  and  Vj  as 

K-  1 

A  =  2<TjU;v/.  (B6) 

;=0 
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Using  this  expression,  a  reduced  rank  matrix  B  can  be  formed  simply  by  setting  the 
appropriate  number  of  scalar  coefficients  to  zero. 

To  simplify  the  problem,  it  is  assumed  that  the  values  of  are  ordered  so  that 
for  a  rank  reduction  of  1,  ck-i  =  0,  for  a  rank  reduction  of  2,  cTjf-i  =  (^k-2  =  0,  etc. 
Examining  the  case  where  the  rank  is  to  be  reduced  by  1,  the  reduced  rank  matrix  may  be 
described  oy, 


K-  2 

B  =  (B7) 

;=0 

Using  equation  (B6)  for  the  definition  of  A,  and  equation  (B7)  for  the  definition 
of  B,  then  the  error  matrix  defined  by  equation  (Bl)  becomes 

K-i  K-2 

E  =  5^<ryuyv/  -  5^a;U;V;”.  (B8) 

;=0  ,=0 

and  simplifying 

,  E  =  ^rjr-iUjr-iUjr-i  (B9) 

If  Ui  is  an  element  of  the  vector  u/r-i  and  v,  is  an  element  of  the  vector  v/m,  then  the 
elements  of  the  matrix  E  are  given  by, 

fit;  =  <rK-lUiV, 

and  the  error  power  is  given  by, 

m-  1  ti  -  1 

{<TK-lUiVj){aK-\UiVj).  (BIO) 

i  »  0  ;■  *  0 

This  expression  can  be  simplified  by  first  rearranging  the  coefficients, 

m  n 

«  *1  ;*l 

H  H 

and  then  recalling  that  v/r-ivir-i  =  Uir-iUjM  =  1,  therefore 

=  ck-\-  (B12) 

2 

Using  this  last  result,  the  error  power  will  be  minimized  if  a  k-i  is  minimized.  In 
terms  of  the  choice  of  the  basis  vector  Vir-i,  and  utilizing  the  relationsUp  in  equation  (B5), 
equation  (B12)  rewritten  as, 


.2- 


K  H  H 

=  o-l-i  =  (<TAr-iii;n)  (ajr.iii^i)  =  Vjr-1  A  Avjr-i. 


(B13) 
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Inspection  of  equation  (B13)  suggests  an  alternate  description  of  the  problem.  That  is, 
choose  a  direction  v^-i  in  the  vector  space  defined  by  the  rows  of  A  (other  than  the  null 
vector)  for  which  the  sum  of  the  power  of  the  row  vector  components  (7j/r-i  in  equation 
(B4))  in  that  direction  is  minimized. 

> 

Minimizing  ^ith  respect  to  each  of  the  elements  of  v/f-i  with  the  constraint 
that  V/r-i^Vjr-i  =  1  is  simplified  using  the  Lagrange  multiplier  technique  [5-8].  Defining  a 
new  function,  *■ 


F(vk-i)  =  vjf-i"AVvir.i  +  Air-i(l  -  Vir-iV-i),  (B14) 

this  incorporates  both  equation  (B13)  and  the  unit  length  constraint  for  Vjf.i.  The  solution 
is  derived  by  minimizing  this  equation  with  r^pect  to  each  of  the  elements  of  v«:.i.  The 
result  is. 


2 A  Avjf-i  -  2A|f.iv/f.i  =  0.  (B15) 

Rearranging  gives, 

AVvjf-i  =  A/r-iV/f-i,  (B16) 

•  H 

which  is  an  eigenvector  expression  for  A  A. 

Returning  once  again  to  the  error  power  expression,  equation  (B13),  and 
incorporating  the  latest  result,  then 

Crain  =  v/f-iV^AvA--!  =  v/r-i^A/f-iVjr-i  =  A/r-i-  (B17) 

Therefore  the  minimum  error  power  Cmin  results  if  the  eigenvector  associated  with  the 
smallest  eigenvalue  of  the  matrix  A^A  is  used  for  the  basis  vector  Vfc-t. 

Extending  this  analysis  to  the  case  where  the  rank  of  matrix  A  is  reduced  by  any 
number,  it  is  apparent  that  the  eigenvectors  of  the  matrix  A”A  should  be  ordered  so  that 
the  corresponding  eigenvalues  are  arranged  in  decreasing  order,  that  is,  Ao  >  Ai  >  A2  >  ...  > 
Ajf-i.  The  optimum  reduced  rank  matrix  can  then  be  computed  by  using  equation  (B7)  and 
setting  the  limit  of  the  summation  equal  to  the  rank  of  the  reduced  matrix. 

From  the  preceding  analysis  it  is  apparent  that  the  basis  vectors, 

▼0,  ▼!,  vj,  Vir-i,  can  be  computed  by  performing  an  eigenanalysis  of  the  matrix  A”A. 

The  corresponing  vectors  uo,  ui,  a2, ...,  U/f-i,  and  scalar  values  ao,  ai,  a2,  aK-i  can  then 
be  computed  using  the  relationship  expressed  by  equation  (B5).  Alternatively,  noting 
vector  set  Uo,  ui,  ...,  u,^i  is  orthonormal  (proof;  Equation  (B5)  can  be  modified  as, 
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If  j  then 


and  if  i  =  j  then 


(cT.U,)“(a,U;)  =  (Av.)\Avj) 

0-,(r;U,”u;  =  V,  “(A”AT;) 

H  H 

U;  =  V,  (AjV  j) 


H 

U,  U; 


Ju_ 

<7,  <Tj 


V, 


H 


H  A  .  H  A .  « 

H  Ai  H  Ai/,\  1 

"5? =  (*)  =  *> 


(B19) 


(B20) 


(B21) 


where  the  relationship  A,  =  (T?  follows  from  equations  (B12)  and  (B17).),  and  modifying 
equation  (B5)  in  the  following  manner: 


(Av;)“=  (a;iij)" 

H  .  H  H 

V  J  A  =  0  jU; 

H  .H  ,  H  V  H 

Vj  A  U;  =  <r;(V;  ▼>;  U; 

=  (B22) 

By  symmetry  then,  the  vector  set  uo,  ui,  U2,  ...,  U/r-i  is  an  orthonormal  basis  vector  set  for 
the  columns  of  matrix  A  [6-3].  Following  the  same  analysis  as  previously  then  U;  can  be 
shown  to  be  an  eigenvector  of  the  matrix  AA"  with  the  same  non-zero  eigenvalue  A^  as 
the  matrix  A”A. 

Given  that  both  vector  sets  uo,  ui,  U2, u/r-i  and  vo,  vi,  V2,  •••,  vx-i  are 
orthonormal  vector  sets,  and  the  values  (Tq,  <j\,  (Tk-\  are  positive  real  numbers,  then 

equation  (B6)  defines  the  singular  value  decomposition  of  matrix  A.  In  this  case  the 
vectors  uo,  ui,  U2,  UiM  are  called  the  left  singidar  vectors,  the  vectors  vo,  vi,  ...,  vjf-i 
are  called  the  right  singular  vectors,  and  the  positive  real  scalar  values  ao,  au  (^2,  ■  <tk-i 
are  called  the  singular  values.  Using  this  fact,  the  singular  vectors  and  singular  values 
may  be  calculated  directly  using  singular  value  decomposition  techniques.  Again  for  the 
purposes  of  optimum  reduced  rank  matrices,  the  ordering  of  the  vectors  will  be  such  that 
for  the  corresponding  singular  values  <7o  >  <ri  >  <7^2  >  ...  >  crg-i  (since  <7^^  =  A;  as  established 
by  equations  (Bl2)  and  (B17}). 
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